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Various theoretical methods address transport effects in quantum dots beyond single-electron 
tunneling, while accounting for the strong interactions in such systems. In this paper we report a 
detailed comparison between three prominent approaches to quantum transport: the fourth order 
Bloch-Redfield quantum master equation (BR), the real-time diagrammatic technique (RT) and the 
scattering rate approach based on the T-matrix (TM). Central to the BR and RT is the generalized 
master equation for the reduced density matrix. We demonstrate the exact equivalence of these two 
techniques. By accounting for coherences (non-diagonal elements of the density matrix) between 
non-secular states, we show how contributions to the transport kernels can be grouped in a physically 
meaningful way. This not only significantly reduces the numerical cost of evaluating the kernels, 
but also yields expressions similar to those obtained in the TM approach, allowing for a detailed 
comparison. However, in the TM approach an ad-hoc regularization procedure is required to cure 
spurious divergences in the expressions for the transition rates in the stationary (zero-frequency) 
limit. We show that these problems derive from incomplete cancellation of reducible contributions 
and do not occur in the BR and RT techniques, resulting in well-behaved expressions in the latter 
two cases. Additionally, we show that a standard regularization procedure of the TM rates employed 
in the literature does not correctly reproduce the BR and RT expressions. All the results apply 
to general quantum dot models and we present explicit rules for the simplified calculation of the 
zero-frequency kernels. Although we focus on fourth order perturbation theory only, the results and 
implications generalize to higher orders. We illustrate our findings for the single impurity Anderson 
model with finite Coulomb interaction in a magnetic field. 

PACS numbers: 73.23.Hk, 73.63.-b, 73.63. Kv, 73.40.Gk 



I. INTRODUCTION 

The experimental progress in fabrication of ultrasmall 
electrical devices 1 - - — has made quantum dots one of the 
standard components in fundamental research and appli- 
cation oriented nanostructurcs. Whereas high-resolution 
transport measurements in the low temperature regime 
have reached a high degree of sophistication and reveal 
data dominated by complex many-body phenomena^—, 
theoretical methods are still struggling to describe these, 
mainly due to competing influences of strong local inter- 
actions and quantum fluctuations^—. A common setup 
for transport studies is drawn in Fig.[TJa): the number 
of electrons on the device is controlled capacitivcly via a 
gate voltage V g , a difference in the electro-chemical po- 
tentials of the leads is created by a bias voltage V\y. The 
measured quantity is the current I or the differential con- 
ductance d//dVb of the whole circuit, which is usually 
represented in stability diagrams, where the changes in 
current vs. V g and Vt, are color-coded [Fig. QJb)]. The 
observed phenomena strongly depend on the strength of 



the coupling of the nanodevice to the electronic reser- 
voirs. In the limit of extremely weak coupling, current 
at low bias is completely blocked in wide ranges of the 
gate voltage, showing up as so-called Coulomb diamonds 
[e.g. Fig. [lib), central region]. Outside these regions 
of Coulomb blockade only single electrons can be trans- 
ferred sequentially onto or out of the do t 34 ' 35 , called 
single-electron tunneling (SET). For simple systems (e.g. 
without orbital degeneracies) in this regime rate equa- 
tions 3 ^ are the standard technique to calculate the oc- 
cupations of the dot states, the current and other trans- 
port quantities 3 -!. The transition rates are calculated by 
Fermi's Golden Rule, i.e., leading order perturbation the- 
ory in the tunneling. For more complex quantum dots 
with degenerate orbitals 3 -^— and/or non-collinear mag- 
netic electrodes^—, coherences, i.e., non-diagonal den- 
sity matrix elements, give crucial contributions to the 
transport quantities and cannot be neglected. These are 
typical situations in molecular electronics and spintron- 

iC9« 

Since the transparency of the contacts is a matter of the 
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FIG. 1: (a) Typical transport measurement circuit setup. The 
source-drain bias voltage Vt, drives a current / through the 
quantum dot. Applying a voltage V s to the gate electrode (in 
this case a backgate beneath an insulating substrate) shifts 
the energy required to add an additional electron to the dot 
by — eaVg, where — e is the electron charge and a the gate 
coupling, (b) The measured differential conductance is usu- 
ally presented in the form of a stability diagram, i.e., in a gate 
voltage - bias voltage plane. Here we show the characteris- 
tics for a single interacting, spin-split level, i.e., an Anderson 
model in magnetic field, which we return to in Sec. [V] and 
Sec. |Vl] (see Fig. H for details). 



material choice as well as of fortune, on the way to low 
ohmic contacts, intermediate coupling strengths are of- 
ten observed, allowing for coherent tunneling of multi- 
ple electrons^. Also, there is the possibility to design 
structures with tunable tunnel barriers, such that dif- 
ferent coupling regimes can be systematically accessed, 
allowing for more detailed spectroscopic information to 
be extracted 8 -. Therefore, electron transport theory must 
go beyond lowest order perturbation theory in the tun- 
neling, while including many charge states, their com- 
plex excitations and their quantum coherence. In recent 
years, several advanced approaches that address higher- 
order effects have been developed based on iterative real- 
time path-integral methods^, scattering-states^ com- 
bined with quantum-monte Carlo^ or numerica l 55 i 56 or 
analytical renormalization group method s 22 i 29 . Although 
these new methods are promising, the standard general- 
ized master equation (GME) approach still offers several 
advantages. The GME describes the reduced density ma- 
trix of the quantum dot with transport kernels that are 
calculated perturbativcly. The GME can be derived us- 
ing various methods^: the Nakajima-Zwanzi g 58 ' 59 pro- 
jection operator technique^— , the real-time diagram- 
matic technique2^£i (RT) and the Bloch-Redfield ap- 
proach^— (BR). Evaluating the kernels up to fourth 
order in the tunneling Hamiltonian (next-to-lcading or- 
der), one can account for all processes involving coher- 
ent tunneling of one or two electrons. Thus, systematic 
corrections to SET can be evaluated either analytically 
in simple cases, or, in complex cases, in a numerically 
efficient way. Although the GME is clearly limited to 
moderate values of the tunnel coupling compared to tem- 
perature, it has the further benefit of non-perturbatively 
treating both the interactions on the dot as well as the 
non-equilibrium conditions imposed by the bias voltage. 
It can therefore provide crucial physical insights into 



measurements of non-linear transport through complex 
quantum dots, see e.g. Re fi 68 ' 69 . 

This paper focuses on the fourth order BR and RT 
methods applied to transport and addresses crucial tech- 
nical issues and simplifications relevant for the descrip- 
tion of complex quantum dots. Several important con- 
crete issues have motivated this work: 

(i) There is an ongoing discussion about the validity 
and equivalence of approaches which can become ob- 
scured by the complexity of the expressions involved 
when discussing complex quantum dots. Clearly, the gen- 
eral form of the quantum master equation is well known 
since several decades. Still, a much debated issue is 
the actual task of systematically calculating higher or- 
der corrections to the transport kernels occurring in this 
equation for general complex quantum dot models. This 
paper emphasizes that the BR and RT techniques are 
one-to-one equivalent. In contrast, the scattering rate 
approach based on the generalized Fermi's golden rule 
and T-matrix (TM), as formulated in the literature, dif- 
fers from these two techniques^. Although it also relies 
on a fourth order perturbative calculation, the results do 
not coincide in general for identical models. The reason 
is that the objects calculated perturbatively, i.e. the T- 
matrix and the time evolution kernel, respectively, are 
different objects whose relation needs to be clarified. In 
particular, the divergences occurring in the TM method 
are intrinsic to the method and not to the problem. We 
show that these go back to a lack of cancellation of diver- 
gent, reducible contributions to the transport kernels and 
that the regularizations proposed in the literature can- 
not reproduce the exact fourth order kernel. We quan- 
titatively demonstrate the resulting deviations from the 
correct GME (BR or RT) result for the example of a 
single impurity Anderson model in magnetic field and 
analytically show how the divergent TM expressions are 
automatically regularized in the GME approaches. The 
GME approaches consistently account for all contribu- 
tions to the perturbation expansion of the transport ker- 
nels in a given order. The importance of this was recently 
highlighted for the well studied non-equilibrium Ander- 
son model, which was found to exhibit a previously un- 
noticed resonance due to coherent tunneling of electron 
pairs^I. 

(ii) The importance of non-diagonal elements in lowest 
order calculations involving degenerate states has long 
been recognized ( "secular contributions" ) , and continues 
to attract attention in the context of transport. Only re- 
cently, the importance of non-secular terms (coherences 
between non-degenerate states) was found to be crucial 2 ^ 
for fourth order tunnel effects. We generalize the dis- 
cussion in Refi2£ and show how these non-secular cor- 
rections can efficiently be included into effective fourth 
order transport kernels through "reducible" diagrams. 

(iii) Explicit expressions for the fourth order transport 
kernels for a very general class of quantum dots were de- 
rived in Refill. However, the numerical cost of evaluat- 
ing these expressions limits their applicability to systems 
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where a relatively small number of many-body excita- 
tions (< 100) has to be accounted for. Here we show 
how contributions to the effective kernels can be grouped, 
making generally valid cancellations explicit and result- 
ing in fewer and simpler terms in the perturbation ex- 
pansion. From direct comparison between numerical im- 
plementations of the expressions in Refill and of our new 
"grouped" expressions, we find the latter to be between 
10 and 20 times faster, without introducing any additional 
approximation. This allows the treatment of more com- 
plicated and realistic quantum dot models. The direct 
gain due to the reduction of the number of expressions 
amounts to a decrease of the computation time by a fac- 
tor of 4. However, the grouping structure can be ex- 
ploited further to make the numerical implementation 
more efficient, leading to the additional speed-up. More- 
over, the grouping gives a basis upon which an explicit 
connection to TM expressions can be revealed. 

As will become clear in the course of this paper, our 
newly found grouping intimately connects and enlightens 
the above three issues, which warrants our systematic 
and extended discussion. 

The structure of the paper is as follows. In Sec. |TT] 
we discuss the model Hamiltonian of the setup Fig. Ufa) 
and some pertinent notation. We then introduce the re- 
duced density matrix (RDM) describing the quantum dot 
as part of the whole system and the generalized master 
(or kinetic) equation (GME) which describes its time- 
evolution. We summarize its general properties and the 
common ground of the discussed approaches. We then 
turn to the derivation of the generalized master equa- 
tion using the BR and RT technique. The crucial role 
played by time-ordering, irreducibility of contributions 
and analytic properties (lack of spurious divergences) is 
emphasized. The derivations are given as compactly as 
possible, because there exists a broad formal study on 
different master-equation approaches by Timrn^i, who 
also showed their equivalence. In contrast to his work, 
our comparison continues in Sec. IIIII on a more explicit 
level with a mapping between terms arising from the BR 
and RT diagrams. In Sec. IIVI we demonstrate how the 
non-secular contributions of coherences lead to important 
corrections in the fourth order transport rates. Based on 
this, Sec. IVl introduces a grouping of contributions to the 
transport kernels, yielding significant simplifications due 
to partial cancellations, which is followed by an analysis 
of how the groups of diagrams contribute to fourth or- 
der physical transport processes (cotunneling, pair tun- 
neling, and level renormalization and broadening). In 
Sec. I VII the derivation of the TM is reformulated. The 
long-standing problem of a precise comparison with the 
RT technique in the context of transport theory and the 
origin of divergences in the TM approach is solved for the 
general case. The theoretical discussion in the last two 
sections, [V] and IVI1 is illustrated by the tangible appli- 
cation to a single impurity Anderson model in magnetic 
field. 



II. MODEL AND GENERALIZED MASTER 
EQUATION 

The standard model for a quantum dot system coupled 
to contacts reads 



H + H-t 



R ■ 



The Hamiltonian 



H R = 



\aq 



w)4 



aqC-lcrq 



(1) 



(2) 



l—s,d 



models the reservoirs, i.e., the source and the drain con- 
tact. The operator c[ (q ct «j) creates (annihilates) an 
electron in a state q with kinetic energy ei aq in the 
source (I = s) or drain (l = d) contact, where a de- 
notes the spin projection. The bias voltage shifts the 
electro-chemical potentials of the source and drain leads 
such that fi s — fid = eVb, where — e is the electron charge. 
The coupling between the quantum dot and the leads is 
described by the tunnel Hamiltonian 



; aqm 



d(jm C lcrq 



(3) 



where d\ m (d am ) creates (annihilates) an electron in the 
single particle state m on the dot. The single-particle 
amplitude ti mq for tunneling from an orbital state q in 
lead I to an orbital state m on the dot is assumed to be 
independent of spin. Finally, the dot is described by the 
Hamiltonian 



H = J2E a \a)(a\, 



(4) 



where \a) is a many-body eigenstate of the dot with 
energy E a . The precise dependence of these energies 
on the applied voltages arising from capacitive effects 
(see e.g., Re fi 35 ' 70 ) is irrelevant for the following discus- 
sion. Typically, the gate voltage dependence is linear, 
E a oc — eV g N a , where N a is the number of electrons for 
state a. 

The diagonalized many-body Hamiltonian Q together 
with the tunnel matrix elements (TMEs) Tj^ q (a,a') be- 
tween all the many-body eigenstatcs a, a', 

T il q ( a , a ') : = £ iim «( a l d LJ a ')' ( 5a ) 

m 

T,- q M := [T+ q (a\a)}\ (5b) 

form the crucial input to the GME transport theory, 
which thereby incorporates local interaction effects non- 
perturbatively. Here Tj^ q (a,a') is only nonzero if the 
state a differs from a' in its electron number by N a — 
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N a > = ±1. Typically, the TMEs can be assumed inde- 
pendent of q, but this is not a prerequisite for the results 
of this work. 

We focus on the regime of weak tunnel coupling and 
thus split H tot in a free part Hq = H + H R and a pertur- 
bation Hj. The condition for weak coupling is that the 
broadening hTi induced by tunneling processes from and 
to lead I is small compared to the thermal energy, i.e., 
hYi <C fcflT, where K and Ub are respectively the Planck 
and Boltzmann constants. The broadening is defined by 

I 1 / = ^ ^2 l^mg| 2 S{ei aq - w), (6) 

where for the contexts in which we need it here (i.e., as a 
measure of magnitude), both the orbital (to) dependence 
and the frequency (u>) dependence arc neglected. Notice 
that the analytical expressions we derive in this work are 
a priori not subject to these restrictions. For simplicity, 
we will further denote contributions of the order of Ti 
in general by T. In this paper we will go beyond lowest 
order in T, allowing the regime of intermediate coupling 
to be addressed. 

Throughout the paper, we will make use of the Liou- 
ville (supcropcrator) notation, in addition to standard 
Hamiltonian operator expressions. The former is merely 
an efficient bookkeeping tool on a general level, whereas 
the latter may be more useful when evaluating matrix 
elements. For example, for the dot Hamiltonian the ab- 
breviation 

CB:=l[H,B] (7) 

defines the action of the Liouville supcropcrator in the 
Schrodinger picture on an arbitrary operator B. It gener- 
ates the time-evolution through e lCt B = e^ m Be~^ m = 
B(t) (Bakcr-Campbcll-Hausdorff formula). Analogous 
expressions hold for the other Hamiltonians, and in par- 
ticular we will need 

rt(t)B I (t') = ±[H$(t),B I (t')}, (8) 

where B 1 \t') = e i(H+H R )t' Be ~i(H+H R )t' = e i(C+C K )t B 
is an operator in the interaction picture. Notice that the 
time evolution of the superoperator in the interaction 
picture is thus 

£{(t) = e l{c+CR)t C T e -^ c+c * ]t . (9) 

A. Generalized master equation and steady state 

The object of interest here is the reduced density ma- 
trixli (RDM) 

p(t)=Tr R {p tot (t)}. (10) 



It describes the state of the quantum dot incorporating 
the presence of the leads, which are traced out of the 
total density matrix p t ot, as prescribed by Tr r. Once p(t) 
is known, the expectation value of any observable can be 
calculated, as discussed below. Before the interaction Hj 
is switched on at time t = to, the total density matrix 
Ptot is the direct product of the (arbitrary) initial state 
p(to) of the quantum dot and the equilibrium state p R of 
the leads, 



e -H R /T 



with Z R = TrR/jR. After the interaction Hj is switched 
on, i.e. for times t > to, correlations, which are of the or- 
der of the tunnel coupling^, build up between leads and 
quantum dot, causing p toi to deviate from the factorized 
form: 



= p(t)p R (t) + 9(t-t Q )0(£ T ). (12) 



We emphasize that it is crucial to include in a kinetic 
equation for p(t), the correlations O(Cj) ~ 0{Hj) be- 
tween leads and quantum dot consistently beyond linear 
order in Hj, if one is interested in going beyond low- 
est order. As we will see in Sees. Ill Bl and III C[ the RT 
approach incorporates them automatically by directly in- 
tegrating out the leads for times t > to, while for the BR 
one explicitly solves for the deviation from the factorized 
state. Both the BR and the RT technique lead to the gen- 
eralized quantum master equation (or kinetic equation), 
describing the time evolution of the RDM 



p(t) = -i£p(t) + [dr K(t - r)p(r). (13) 

Jt 



Here, the first term accounts for the time evolution due 
to the local dynamics of the quantum dot. In the second 
time non-local term, the time evolution kernel K.{t — t) 
is a supcropcrator acting on the density operator. Con- 
voluted in time with p(r), it gives that part of the time 
evolution which is generated by the tunneling. We note 
that this form of the GME is dictated by the linearity of 
the Liouville equation and the partial trace operation. 

The kernel to fourth order formally reads 
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K(t - t) P (t) = — Tr r {£Te- i(£+£R)(t - r) £Tp(T)p R } 



(14) 



J dr 2 dn [Tr R {£ T e- i;(£+£R)(t - T2) /:Te-' i(£+£R)(T2 - Tl) /:Te- i;(£+£R)(Tl - T) /:Tp(T)pR} 

*>T 2 >Ti>T 

—TrR [Cje-^+^-^Cje-^+^-pR Tr r {e^+^^Te-^+^^-^xKr)/*} } 

I 



Below we will show how BR and RT (as well as the 
Nakajima-Zwanzig projection operator approach, see 
App.[A}, when consistently applied, lead to this result. 
The central topic of this paper is the explicit evaluation 
of Eq. (|14p and the significant simplifications which can 
be achieved in the steady state limit lim^oo p(t) = p. In 
this limit, Eq. (fT3]) becomes 



(15) 



lim pit) = = — iCp + Kp, 

t—>oo 



where K = K(z = iO) and iO denotes an imaginary in- 
finitesimal, and 



/>oo 

K(z)= / dre lZT IC(T) 
Jo 



(16) 



is the Laplace transform of the time evolution kernel. 
Taking matrix elements with respect to the many-body 
cigenstates of the dot Hamiltonian, H, we obtain from 
Eq. ([15)) a set of linear coupled equations for all states 
a, a' of the RDM: 

= -t <W (E a - E a ,) p aa , + ^ K$ ' Pa a> ■ (17) 



The expectation value of any non-local observable can 
be expressed in a form similar to Eq. (|13p . In particular, 
we can write the particle current flowing out of lead I (i.e. 
the number of electrons leaving lead I per unit time) as 



Il(t) = (Ii(t))=Tr / drAC/,(t-r)p(r) 



(21) 



where ICi l (t — r) is the kernel associated with the current 
operator 



h = 



- t [Hti,Ni] 
n 



= -zX! tlmqdlmClvq + h - C ' 



(22) 



aqra 



with Ni = J2aq c \aq c i°q being the number operator in 
lead I. 

Taking the steady state limit of Eq. (|2"Tj). the DC current 
is given by the zero- frequency component Kj t := Kj t (z = 
iO) of the Laplace transform of the current kernel, traced 
over with the stationary density matrix p: 



I t =Tr {K h p} = Y d {K I 



\ aa n , 
Ibb Paa' 



(23) 



Here, the matrix elements of K (or any other supcropcr- 
ator) are defined by 



Kw -=(b\[K\ 



')(< 



16') 



(18) 



where we use square brackets to make clear that the ker- 
nel superoperator must first act on |a)(a'|, and then the 
matrix elements of the resulting operator are taken. Each 
diagonal element of the RDM reflects the probability of 
finding the system in a certain state. Thus, the normal- 
ization condition 



(19) 



must be fulfilled. The restriction (jl9|) allows the system 
of linear equations obtained from Eq. (jTTJ) to be solved, 
since without it they are under-determined due to the 
sum-rule 



E 



rs-aa 
^bb 



Ma,a'. 



(20) 



Physically, this guarantees that gain and loss of proba- 
bility arc balanced in the stationary state. 



The current kernel can be calculated by simple modifi- 
cation of the time-evolution kernel as discussed below in 
subsections III Bl and III CI explicitly. We will now address 
the derivation of Eq. (|13J) and of its kernel (TT4"|) up to 
fourth order in the tunnel coupling. We focus here on 
two approaches, an iterative procedure in the time do- 
main^r— , referred to as Bloch-Redfield approach (BR), 
and the real-time diagrammatic techniqu o 25 ' 63 ' 64 (RT). 
The projection operator technique of Nakajima^ and 
Zwanzig^, which has been explained and used in many 
works££ ~ 62 i 72 , is closely related and equivalent to the BR 
approach and is discussed for completeness in App. |XJ 
The derivation of the kinetic equation requires no other 
ingredient than the Liouville equation for the total den- 
sity matrix p to ^' 



p I tot (t) = -t£ I T (t)pUt). 



(24) 



For the purposes of this paper, it is most convenient to 
work in the time-domain and use the interaction picture. 
In addition, we make use of the property of the particular 
bi-linear coupling of Eq. considered here, that the 
lead-average of an odd number of interactions vanishes 
due to the odd number of lead field operators in Hj. 
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B. Bloch-Redfield approach 

The Bloch-Redfield approach^— is usually favored to 
derive the second order quantum master equation^. Ba- 
sically, one integrates Eq. (|24p and reinserts it back into 
its differential form to get 

pUt) = -iC{(t)pL(to) - [ dr C{{t)C{{T)p{ ot {T). 

(25) 

We now extend this to fourth order— by repeating the it- 
eration steps: we transform Eq. (|25[) to an integral equa- 
tion, 

PtotW = pLih) -if dr £{(T)pUtv) 

J t 

I dndr#(n)#(T)^(T), (26) 

t>r 1 >T>t 

which is once more reinserted into Eq. (|24[) . After inte- 
gration one arrives at 

PUt) = Ptot(io) - i /dr C{(r)p( ot (t ) 

J to 

- J dndrC{(n)C{(r)pUto) 

t>T t >T>ta 

+ i J dndndrjC^C^C^pUr)- (27) 

t>T2>T 1 >T>t 

We reinsert Eq. (|2~T)) back into the Liouville equation (f2"4")l 
and perform the trace over the leads in order to ob- 
tain the RDM. Thereby, terms involving in total an odd 



number of lead operators vanish. Due to the relations 
p( ot (to) = p I (to)pR and with Eq. (TT2"]) we obtain: 

p\t) = - I dr 2 Tr R {4(i)4(r 2 )p I (t ) / 9 R } 

Jto 

+ J dradndrTrp {L^C^jC^C^p 1 (r)p R } 

*>T2>Ti>T>i 

+ 0((C{n (28) 

The second order contribution in Eq. ((28)) contains p 1 (to) 
instead of p 1 (r) and lacks the convoluted form which the 
fourth order term has: thus in the stationary limit the ini- 
tial state p 1 (to) does not seem to drop out. If one naively 
were to neglect this difference and set p 1 (to) ~ p 1 (t) the 
fourth order kernel would contain spurious divergences 
(see Sec. |VI| . Instead, one has to account for the corre- 
lations between dot and reservoirs at times t > to up to 
order (£{) 2 as expressed by Eq. (|2"6"|) . Taking the trace 
over the leads this equation gives: 

/(r 2 ) =/(*o) - J dr.dr Tr R {jC{(n)£{ (t)/(t)Pr} 

T2>T 1 >T>t 

+ 0((C{f). (29) 

This shows also that Eq. (|28|) through p 1 (to) still con- 
tains higher-order terms. Consistently neglecting these 
in Eq. (|28p . we can eliminate the dependence on the ini- 
tial condition from Eq. (f2"5|) and thus arrive at the GME 
in the interaction picture 

p I (t)= f dr /C 7 (^t)/(t), (30) 

Jto 

with the time-evolution kernel defined by 



JC'(t - r)/(r) = —TrR {4(t)£((r) /9R p / (r)} + 

J dr 2 dn [Tr R {£((i)£((r 2 )4(Ti)4(r) /9R /(T)} -Tr R {£|(t)£|(r 2 ) p R Tr R {£f (n)£f (r)p R /(r)} }] . (31) 

t>T 2 >T 1 >T 

I 



Transforming the RDM to the Schrodingcr picture by 



p I (t) = e^+^ t (p(t)+iCp(t)) 1 (32) 



and with the Liouville operators according to Eq. ([9]), 
we obtain the generalized master equation (|13[) . and 
arrive at the expression (|14[) for the kernel. 



The current, Eq. (f2~Tj) . is given by 

Il(t)= Tr tot {ll(t)pUt)} (33a) 
= Tr / dt' K} h (t - t^p 1 ^'), (33b) 

Jto 

where the current kernel in the interaction picture is 
given by Eq. (f3Tj) with £{(t) replaced by h(t). In de- 
riving this, as for the density matrix, one must take 
care to keep the time-ordered structure: since the cur- 
rent operator h(t) is, as Hji, linear in the lead opera- 



7 



tors [cf. Eq. (|22[) ]. we obtain Eq. (|33bp , correct up to 
fourth order, by inserting the third order iteration for 
Ptot (Eq. p7]) l into Eq. (|33a|) . Under the trace, the first 
and third contribution are zero since they contain an odd 
number of interactions. Next, as before, in the second 
contribution of Eq. (f27j) . p(£o) has to be eliminated us- 
ing Eq. (|29| , thereby generating a fourth order correction 
term. Finally, in the fourth term one must consistently 
keep ptotij) ~ p(t)pr, i.e. only here one can neglect the 
deviation from the factorized form. 



C. Real-time diagrammatic technique 

The real-time approach has been discussed on a gen- 
eral level in many work o 25 ' 63 ' 64 . Therefore the aim of 
this section is to recall how one efficiently arrives at 
the kinetic equation and its kernel by exploiting Wick's 
theorem from the outset. We start from the Liouville 
equation (|24|) for the full system in the interaction pic- 
ture and formally integrate it: 

P Ut)=Te- i ^ drC ^pUto), (34) 

where T is the time-ordering superoperator. Using 
/Otot(*o) = Ptot(to) = Pr(*o)/' I (*o), and defining the su- 
peroperator 

IT'&to) = Tr R {Te- lJ l dTC ^ T} PR (to)} , (35) 

the time-evolution of the reduced density matrix can for- 
mally be written as: 

/(i) = 7T I (t,t )p I (t ). (36) 

Expanding the time-ordered exponential superopera- 
tor, the trace in Eq. (|35[) can be explicitly evaluated term 
by term by Wick's theorem: the trace over each string of 
reservoir field operators becomes a product of pair con- 
tractions, indicated in the following by contraction lines. 
For our purposes here, one can simply formally consider 
the Liouvillians to be contracted (meaning their reser- 
voir part), see Eq. (|38|) below. We can then decompose 
7T (t, to) into a reducible and an irreducible part, de- 
pending on whether or not the contractions separate into 
disconnected blocks. Collecting all irreducible parts into 
the kernel K, 1 {t - t') one obtains in the standard way a 
Dyson equation: 

7T J (Mo) = 1 + J dr 2 dn K\n - t^TT 1 (r u t Q ). (37) 

*>T2>Ti>*o 

It relates the full propagator 7T 7 (Mo) to the free propa- 
gator, which equals unity in the interaction picture, and 
to the irreducible kernel K} . Applying the Dyson equa- 
tion to p(to) and taking the time derivative, one arrives at 
the kinetic equation in the interaction picture, Eq. (|30[) . 
Transformed back to the Schrodinger picture we obtain 



the kinetic equation (|T3l) . 

We have now obtained the kernel JC 1 (t — t) formally as 
the sum of irreducible contributions to the time-evolution 
superoperators of different orders in the tunneling, which 
are written down to fourth order: 

lC\t -r) = -L{{t)L\{r) + J dr 2 dn 

t>T2>Tl>T 

- I l=i I I i=i I . 

4(i)£((r 2 )^(r 1 )^(r) + ^(t)4(r 2 )^(r 1 )^(r) , 

(38) 

where £f(f)#(r) := Tr R {C{ (t)£f (r)p R } . 

The expectation value of the current (or of any op- 
erator) is obtained in a similar fashion. We first intro- 
duce a superoperator = ■}, which is an anti- 
commutator in contrast to the superoperators governing 
the time-evolution. The current is then expressed as 

I l I (t) = Tr{irj l (t,to)p r (t )}, (39) 
where we introduced a current propagator 

TT*. (Mo) = Tr R {T4(M- 4/ 4 d ^V R (i )} , (40) 

which differs from the propagator for the reduced density 
operator only by the current vertex (f ) at the final 

time. Collecting all parts of 7Tj which are irreducibly 
connected to the latter vertex, one readily verifies that 
the remaining irreducible parts at earlier times are those 
contained in the propagator 7T 1 : 

7r^(Mo)= [ dDtdnKfa - n)7T J (TMo). (41) 

f>T2>Tl>to 

The current kernel )Cj (t — r) is given formally as the sum 
of irreducible contributions to the time-evolution super- 
operators of different orders in the tunneling with the 
leftmost interaction vertex replaced by the current su- 
peroperator. Applying this equation to p(tg) and tak- 
ing the trace over the dot, one arrives at the expres- 
sion for the current in terms of the new kernel and 
the reduced density matrix in the interaction picture, 
Eq. (|33b[) . Transformed back to the Schrodinger picture 
we obtain Eq. (|2"T|) . 

D. Comparison of the approaches 

For the comparison between the BR and RT ap- 
proaches, it is most useful to consult Eq. (f3~Tj) and 
Eq. (|38|) . The second order terms, contained in both 
equations in the first line, obviously match. The equiv- 
alence of the fourth-order terms is more indirect: in the 
BR approach, the first term of Eq. ([3"Tj) is evaluated using 
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Wicks' theorem by building all possible contractions, in- 
cluding the reducible ones (contraction between vertices 
at times t and T2 as well as T\ and r). The latter are then 
canceled by the second term. Precisely the same happens 
in the projection operator approach [cf. Eq. (|A6[) ] . The 
above conclusions hold in fact for any order of pertur- 
bation theory as shown inSl. In contrast, the RT ap- 
proach avoids the inclusion and subsequent cancellation 
of reducible parts which rapidly grow in number with the 
perturbation order. 

We emphasize that there is one unique correct fourth 
order (time non-local) generalized master equation, in 
which the kernel includes all fourth order contributions, 
but no higher order ones and which does not divergence 
in the stationary (zero-frequency) limit. This master 
equation can be derived using either the BR or RT (or 
NZ) approaches and there is no need to distinguish be- 
tween these in the following discussion. 

After this comparison on a formal level, we will con- 
tinue in Sec. IIIII with a comparison on the level of the 
individual contributions to the time-evolution kernel. 



III. DIAGRAMMATIC REPRESENTATION 
AND MAPPING BETWEEN BR AND RT 

We now address the task of calculating all elements 
Kgff (iO) of the time evolution kernel in the stationary 
GME, Eq. (|17p. For our purposes, it will turn out to be 
advantageous to first work in the time-domain, i.e., to 
calculate K%fr (t — r), which we decompose into contribu- 
tions of successive non- vanishing even orders n = 2, 4, ... 
in the tunnel coupling: 

={K.W)tf + + ■■• . 

The section has a twofold aim. We first introduce the 
diagrammatic representation for the time evolution ker- 
nel and show how each BR contribution, obtained from 
Eq. (|31[) . translates into a corresponding diagram in the 
RT approach based on Eq. (|38)l . Apart from being of 
technical interest, this is of key importance for the dis- 
cussion in Sec. IIV1 where the fourth order kernel elements 
incorporating corrections to the diagonal elements of the 
density matrix due to the non-diagonal elements are in- 
troduced. Secondly, we discuss the time-dependent part 
of the kernel contributions in the Schrodinger picture and 
its zero-frequency Laplace transform, on which the sim- 
plified calculation of the effective fourth order kernel in 
Sec. IVl relies. 

In the conventional RT approach one starts by consid- 
ering super matrix elements [cf. Eq. (fl"8)) ] of the kernel^ 
and one introduces a diagrammatic representation for the 
order n = 2, 4, . . . parts of the kernel: 



(/CW)£(f-T) 



t ■ 

b 4, 



> 



b' 



The diagram represents operators which act from the left 
and right on the dot density operator p(r), inducing an 
irreducible time-evolution of a pair of initial states a, a' 
[associated with p(r)] to corresponding final states b, b' 
[belonging to p(t)]. Time thus increases from right to 
left. The diagram can be considered as a Keldysh con- 
tour, i.e., running from a — > 6, then continuing backward 
from b' a', as indicated by the directed line on the 
upper and lower part. On the upper (lower) part of the 
contour the time-ordering agrees with (is opposite to) the 
contour direction, indicating that operators acting from 
the right on the density operator p(r) come in inverted 
order concerning time. This distinction is important for 
the diagram rules. The shaded area indicates the sum of 
all contributions involving the product of n tunnel oper- 
ators Ht, starting at time r and ending at time t (this 
thus yields a product of ^ broadening elements, which 
we indicate with T"/ 2 ). Starting from the RT expression 
for the time evolution kernel, Eq. (|38|) . simple rules are 
derived, given in App. IB 2| from which one can directly 
read off the analytical expression for the zero frequency 
Laplace transform of each diagram^ 
Hence, the diagram contains not only the information 
about the contribution to the kernel K, in the interaction 
picture, but also to its Laplace transform K . To make a 
distinction, we will use the convention that we mean con- 
tribution to K, whenever we explicitly write down times 
labels in the diagram. Otherwise it stands for the Laplace 
transformed expression, i.e. the contribution to K (this 
is the case everywhere in the following except for Fig. [2]). 

In the BR approach one has to expand the superop- 
crator expression Eq. (|31[) in commutators with dot and 
electrode operators. One then applies Wick theorem to 
integrate out the electrodes, resulting in cancellations of 
terms. Finally super-matrix elements are taken and the 
remaining expressions correspond term by term to the RT 
diagrams. To emphasize the close connection between the 
two approaches we now illustrate this explicitly by calcu- 
lating second and fourth order contributions to the time 
evolution kernel in the BR approach. To this end, we 
first split the tunneling Hamiltonian (|3|) into two parts, 



A 



k ■ 



describing tunneling into (p — +) and out of (p = — ) the 
dot: 



Cj- := (R- X 4m(n)) (t lmq c laq (n)) , (42a) 
A- k =: C+Dl := (t* lmq cj^fa)) (V^r,)) . (42b) 



Here the index k numbers the time argument, with t and 
r corresponding to k — 3 and k = 0, respectively. The 
summations over l,m,o~,q are implicitly understood. We 
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b +- 



(jc(%;(t-r)= |>i>i>t 



V- 



-¥■ a' 



FIG. 2: Time ordering in a diagram associated with a fourth 
order process. Every term arising from the BR approach, 
Eq. (|3ip . can uniquely be translated into a specific diagram. 
This gives a one-to-one mapping between the BR and RT 
approach. While the time-order is crucial for this mapping, 
the resulting diagram can also directly be used to represent 
the Laplace transformed expression (see Fig. [3]). 



insert this into Eq. (j3"Tj) . denoting pi := p 1 (t): 

K,^{t-r)p I {r) = - Tr R [^,[^>,^p R ] 

(43a) 



e{+ ,-} 



£ w (t-T)/(r)= /dr 2 drix 
po,pi,P2,P3e{+,-} t>T2>Tl>r 



Tr R 



-Tr, 



Af,Tr R {[A{\[^,p I p R ]]}p R 



(43b) 



We next expand the multiple commutators and collect 
the fcrmionic operators of the leads, using that they anti- 
commute with the quantum dot operators, C% k D^ k = 
—D^ k Cl k , where pk = —pk- Using the cyclic property of 
the trace and Wick's theorem, the average of the lead 
fermionic operators is expressed as a sum over prod- 
ucts of pair-contractions. This results for Eq. (|43a[) and 
Eq. (|43b[) in the expressions listed on the left sides of 
Fig. [3] For the fourth order, the reducible contractions 
emerging from the last two terms in Eq. (|3 1 [) cancel each 
other. Using the time-ordering of Fig. [21 Fig. [3] further 
gives on the right sides the respective RT diagram for 
each expression. The Hermitian conjugated terms, which 
have been omitted in the figures, correspond to diagrams 
which are vertically mirrored, i.e. all vertices on the up- 
per contour have to be moved to the lower one and vice 
versa. 

The translation from the BR expressions in the inter- 
action picture into the RT diagrams works as follows. 
For each operator D^ k standing on the left (right) of p$, 
draw a vertex on the upper (lower) contour at time Tfc. 
For each contraction {Cf * Cj 3 ) , requiring pi = —pj in 
order not to vanish, draw a contraction line between the 
vertices representing D? 1 and D^ 3 . Notice that the order- 
ing of the C-operators in each contraction is consistently 
incorporated in the diagram: the pair of C operators 
in the contraction have the same time-ordering as the 
corresponding D operators, unless the earlier vertex of 
the two lies on the lower contour (this follows from the 
cyclic permutation under the trace). Similarly, the sign 



of the operator expression, (— 1)™/ 2 +™c+™i ; j g automat- 
ically contained in the diagram through the number of 
contractions n/2 (n=order of perturbation theory), the 
number of crossing contraction lines n c , and the number 
of vertices on the lower contour n\. 

The diagrams listed in Fig. [3] represent expressions 
which are summed over the indices pk = ±. Terms 
with specific values of the p's are represented by dia- 
grams where the contraction lines arc directed by an ar- 
row, pointing towards the vertex corresponding to D + . 
Figure 0] shows an example for the third 4th order dia- 
gram in Fig. [3] From the diagrams it is explicitly clear 
that all contributions which where not canceled are irre- 
ducible: between the first and the last vertex at times t, 
respectively r, there is no time point at which the dia- 
gram could be separated into two parts without cutting 
a contraction line. 

To obtain the stationary kernel Eq. (f31j) required 
in Eq. (fT5")) we first need to transform back to the 
Schrodingcr picture [cf. Eq. |9"])] by inserting 



_ c iH R T k (t) e -iH R T k 
— e C l<rq C > 



/(r fe ) = e lHTk p(T k )e 



-iHr k 



(44a) 
(44b) 
(44c) 



For the further developments in this paper only the 
time-dependent part of the resulting expression, and 
its Laplace transform, is of importance. We can only 
factor out this part after taking super-matrix elements 
[cf. Eq. ([T5)l ] of the kernel contributions with respect to 
the energy-eigenstates of the quantum dot and insert 
complete sets of these states between all the vertex op- 
erators D. The resulting expression is represented by a 
diagram labeled with these dot eigenstates on the con- 
tour, as illustrated in Fig. 2] We calculate its Laplace 
transform with respect to t' = t — r (r„_i = t, To = r), 
collecting for each time Tk all energy contributions into 
one exponential with argument —iAj-Tk/h: 



pOO p 71 — 1 

/ dr'e^ Aor ' / dT n _ 2 ---dTi J]' 

" h— 1 



t>T n -2--->Tl>t — T' 



(45) 



Here the additional Ao contains both the Laplace vari- 
able z = iO and the energy difference E a < — E a of the 
initial states on the upper and lower contour. Trans- 
forming variables to the time-differences between vertices 
Tk '■= Tfc-i — Tk decouples the integrals, showing that the 
energies 5^ = y\_ n A; fully determine the time-evolution 
factor and its zero-frequency transform: 

n—l „oo n— 1 , 

(k(z))z ~ n / &fk e ~ i5kfk =Ut' (46) 

k=\ Jo fe=i k 

This is the form of the zero frequency Laplace transform 
of the time-dependent factor only, as obtained from the 
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+ <(% C7f 
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BR expressions 
(interaction picture) 

2nd order 

(CECl) (p' DlDl) + (CECl) (Dg P 5 Dg) + h.c. 



associated diagrammatic representation 
(Laplace transformed) 



+ (C p 2 Cf 



(cec!) pi ded(de'de - (cicl) 

(cEcl) 
(CECl) 
(c?c!) 
(cici) 
(cEcl) 
(cici) 
(cici) 



u °2 



(cici) 
(cici) 

(CECl) 



u °2 



(cici) 
(cici) 



DlDt pEDE'Df 

dedIpEd^'di' 

DlDlplDlDi 
d{ pEDIDIDI 



DlpEDl D\D\ 

d{ pEdide'de 



4th order 

'ct) P lD{D{DiDl 

'ct)DlDt pi DlDt 
'C!')D P 3 D( plDlDi 
'cE)dIDI po -Df -Dj' 
'<??') £>f' p' DlDlDE 
'ct)DlpEDlD{Di 

'ci^dEpEdi'di'di 

'cl)DlpEDf;DiDl+h. C . 





+ 






+ 






+ 


* 


^ » 


+ 






+ 




e 


+ 






+ 












+ 





+ h.c. p 



+ h.c. p 



FIG. 3: On the left we give the contributions to the time evolution kernel K. 1 ^ and /C' 4 ' as they arise from the BR in the 
interaction picture. Here, p = ±, p = —p. The corresponding diagrammatic RT representations, standing directly for the 
Laplace transformed contributions to and K^' , are shown on the right. The minus sign of the first BR contribution is 

incorporated into the diagram. 



diagram rules in App. IB 21 which is the most convenient 
starting point for our analysis. The energy 8k is obtained 
from a diagram by making a vertical cut through the di- 
agram between times and Tfc_i, and then adding to 
/ subtracting from z the energies of all directed lines 
which hit this cut from the left / right. This includes the 
energies associated with contraction lines as well as the 
upper and lower part of the contour. Figure 0] demon- 
strates how this simple rule works for a specific fourth 
order diagram. A crucial point for the rest of the paper 
is that diagrams which differ by breaking time-ordering 
between the contours, but keeping time-ordering within 
each contour, only differ by the arguments of the time- 
dependent exponentials in Eq. (|46|) . i.e., in the Sk- As 
App. IB 21 shows, the products of TMEs and electrode 
distribution functions and the overall phase factor are 
identical. The simplifications discussed below are thus 
independent of these factors and their precise form needs 
not be discussed further. 




= z 



c,d 



, b' 




S, -+E.-E +co' 

I da 

d n *E-E 

a a 
a 



5 2 "■* E b ,-E a -co+co' 



S 7 ™+ E,,-E -co 

3 D C 



FIG. 4: Example of one possible fourth order contribution 
to the kernel element K^y . Here, w and u/ are the energies 
which are assigned to the fermion lines. Unless the inter- 
mediate states of a diagram are labelled (like for example in 
Sec. IIVI concerning reducible diagrams), it is implicitly meant 
to contain the sum over any intermediate state. As described 
in the main text, the quantities Si can by read off from cuts 
through the diagram. For our example: So = E a i — E a , 
<5i = i0 + E d — E a + io' , $2 — iO + E b i — E a — ui + ui' , 
S 2 =i0 + E v -E c -u. 



IV. COHERENCES AND NON-SECULAR 
CORRECTIONS 

For many simple quantum dot models, selection rules 
deriving from symmetries prevent the occupations of the 



dot states from coupling to the coherences. Whenever 
two states a and a' of the system differ by some quantum 
number which is conserved in the total system (i.e. in- 
cluding the reservoirs), then their coherence p aa ' does not 
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enter into the calculation of the occupancies since 
for all states. The simplest example for such a quantum 
number is the electron number which is conserved for a 
quantum dot coupled to non-superconducting electrodes. 
The total spin-projection is also conserved for unpolar- 
ized or collinearly polarized electrodes. For non-collinear 
polarizations, however, inclusion of the coherences is cru- 
cial in order to capture spin-precession effect s 45 ' 74 . In a 
similar way orbital degeneracies have been shown to af- 
fect the occupations through the coherences! 42 ' 43 So in 
general, coherences cannot be neglected. Making no spe- 
cific assumptions about the coherences, the only selection 
rule we enforce here is the one due to the conservation of 
total charge. In the above mentioned works, the tunnel- 
ing is treated to lowest order (C^) and only non-diagonal 
elements between degenerate states are kept. The latter 
so-called secular approximation is usually phrased as ne- 
glecting the rapidly oscillating terms corresponding to 
coherences between non-degenerate states^ 

In fourth order, however, an elimination of coherences 
between non-degenerate states in the stationary limit re- 
quires an expansion of the effective kernel for the occupa- 
tions. Such an expansion is consistent in the sense that 
the derived effective kernel includes all contributions up 
to fourth order (while a simple neglect of non-secular co- 
herences introduces serious errors of the order £y)i 27 ' 75 

We start by decomposing the density matrix into a sec- 
ular (energy diagonal) part p s and a non-secular (energy 
non-diagonal) part p n . Here p n contains all matrix ele- 
ments p aa i between states with \E a — E a i\ > e n and p s 
all other elements (including the diagonal components, 
a = a' , corresponding to the populations). The cutoff e„ 
should be chosen large compared to the tunnel broaden- 
ing of the states, e n ^> HY, the precise requirement being 
that it should be large enough that the next-to-leading 
order term in the expansion of Eq. (|47[) below is compa- 
rable to a sixth order term, and can thus be neglected. 
Our aim is to eliminate the non-secular coherences p n 
and include their effect as a correction to the kernel de- 
termining the secular part. To this end we write Eq. (|15[) 
in block matrix form, 



iC r 



K^+K^ \( Pt 



where, see Eq. (fTT)). the free evolution of the system in- 
volving (C)^ = (E a > — E a )S a b6 a 'b> is zero in the ns and 
sn blocks by definition. Solving for p n one obtains 

(47) 

which obviously contains all orders in Y due to the in- 
verse. Since by construction C nn ^> KY we can expand 
Eq. (j47|) in HY/£ nn , finding that the lowest order term 
gives corrections to the secular density matrix of order 
r 2 and is thus all that should be kept in a consistent 
fourth order expansion. Inserting this in the equation for 



the secular part of the density matrix we obtain an ef- 
fective stationary kinetic equation for the secular density 
matrix: 



0=1 -iC 



with the effective fourth order kernel 



K 



-(4) 



Here 



K 



eff 



(4) 



Kit ] + K 



(4) 



nn 



(48) 



(49) 



(50) 



is the correction to the secular density matrix due to 
coherences between non-secular states. This makes ex- 
plicit that when going beyond lowest order, the secular 
approximation is no longer valid and also coherences be- 
tween non-secular states have to be accounted for. This 
was shown ir*2I for the special case where the secular part 
is diagonal. Here we extended the derivation to an arbi- 
trary excitation spectrum, where the secular part may be 
non-diagonal, i.e., the effective equation is not a master- 
equation for occupancies. It should be noted that in this 
case the kernel K$V must be calculated including the 
elements which couple to secular coherences. 

For the developments of the present paper it is useful 
to introduce a diagrammatic representation of the non- 
secular correction . We first note that the inverse of 
C nn is related to a diagram without tunneling lines by the 
diagram rules, see App. IB 21 evaluated at zero frequency 
(z = iO): 



— C 



a -4- 



E a > — E a 



■> a' 



Note that this " free evolution" term is always finite since 
the expansion is only carried out in the non-secular sub- 
space where \E a — E a >\ ^> HY and iO can always be 
dropped. Representing a general second order contribu- 
tion diagrammatically as 



(a-< 2 >) 



hi/ 



b' 



r 



the correction term due to coherences between non- 
secular states is given by 

c 



E 

c,d£n 



1* — - 




r 


r 




— — >* 



where the sum is restricted to states which are pairwisc 
non-secular, i.e., c and d with \E ( i — E c \^> KY. 
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Thus, K^' appears as a sum of all reducible fourth 
order diagrams with non-secular intermediate free prop- 
agating states c, d. The evaluation can be performed by 
using the diagram rules, App. IB 21 as for an irreducible 
fourth order diagram. The effective fourth-order part of 
the kernel, determining the secular part of the station- 
ary density matrix through Eq. (|48[) . can thus be calcu- 
lated in the same way as with only the following 
modifications of the diagram rules: (i) diagrams can be 
irreducible and reducible between the first and last ver- 
tex] (ii) the intermediate states of reducible diagrams are 
restricted to non-secular free propagating intermediate 
states; (iii) only secular initial and secular final states for 
the diagrams are possible. We finally note that when cal- 
culating the current using the current kernel, the contri- 
butions of the non-secular coherences can be eliminated 
in exactly the same way. The only modification required 
is to replace the operator acting at the latest time by the 
corresponding current operator. 

We have thus eliminated the non-secular coherences 
from the transport calculation. This effective diagram- 
matic theory for the secular part of the density matrix 
is the starting point for the evaluation of diagrams in 
groups, rather than single ones, which we turn to now. 
This will result in simplifications of the numerical eval- 
uation of the kernels (Sec. [V}, and allow the relation to 
the TM approach to be established fSec. IVII) . 



V. EFFICIENT DIAGRAM EVALUATION 

We now focus on the explicit evaluation of the effective 
fourth order transport kernel (|48p in the zero-frequency 
limit, 



K. 



(4) 
eff 



K<$ + K 



(4) 



(51) 



according to the modified diagram rules derived above. 
For each super matrix element of this kernel [cf. Eq. (|18[)] 
this involves evaluation of all fourth order diagrams for 
all possible combinations of all intermediate states on 
both propagators. Already for quantum dot models with 
a moderate number of states ~ 10 — 100 this comes at 
a high numerical cost. In this section we demonstrate a 
way to reduce the computational cost drastically without 
introducing any approximation. 

In Fig. [5] we show the 24 diagrams representing the 16 
irreducible and 8 reducible diagrams of K$ and 
respectively. These represent all 192 contributions to the 
effective kernel Eq. (j5"Tj) . since we do not specify the di- 
rection of the two contraction lines (p indices) and only 
include one diagram from each hermitian conjugate pair. 
It is only by combining diagrams from both KsV and 

which is necessary as explained in Sec. IIV1 that a 
structure is revealed upon which our efficient diagram 
evaluation is based. The diagrams are sorted in Fig. [5] in 
three steps. 



A.(0) A.(l) A.(2) 



B.(0) 



B.(l) 



B.(0)(s) 



^ O 



B.(l)(s) 



B.(l)(t) 



z 



-a 



A.(0)(s) 


A.(l)(s) 






5^ 

— ^ — *■ 








A.(l)(t) 


A.(2)(t) 


^ 
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V 







B.(2) 



B.(2)(t) 
tj 



C.(0) C.(l) C.(2) 



C.(0)(s) 


C.(l)(s) 














C.(l)(t) 


C.(2)(t) 


















« / 









FIG. 5: The sixteen irreducible fourth order diagrams, to- 
gether with the eight reducible correction diagrams, can be 
sorted by topology into three diagram classes G = A, B, C. 
Within each class, there are three groups G.(0), G.(l), G.(2), 
labeled by the number of vertices on the upper part of the con- 
tour. In G.(l), the latest (leftmost) vertex distinguishes be- 
tween stand-alone diagrams (s) and triple diagram subgroups 
(t). 
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First, there are three diagram classes, Gs {A, B, C}. 
These distinguish the three topologically different ways 
in which the four vertices can be contracted, considering 
them to lie on the contour (i.e. moving vertices between 
the upper and lower part of the contour via the latest 
time does not alter the topology.) 

Secondly, these classes divide into groups G.(x), x £ 
{0, 1, 2}, based on the number of vertices x on the upper 
part of the contour. (Diagrams with x = {3, 4} need not 
be included in Fig. [5] since they are hermitian conjugates 
to diagrams in the x = {1,0} groups.) The classes are 
thus constructed by forming the one-member group G.(0) 
where all vertices lie on one part of contour (here taken 
to be the lower part) and then successively shifting the 
vertices to the opposite part of the contour via the latest 
time point. 

Thirdly, one distinguishes subgroups by the position 
of the latest vertex, being either on the upper or lower 
part of the contour: the groups G.(l) thus divide into 
a subgroup G.(l)(s) (single) and G.(l)(t) (triple) of one 
and three diagrams respectively, whereas in the groups 
G.(0) and G.(2) there exists only the one- respectively 
three-diagram subgroup, such that G.(0) =G.(0)(s) and 
G.(2) =G.(0)(t). 

A key point is that diagrams within a group G.(x) 
give rise to expressions which in the interaction picture 
only differ by their time-dependent factor, and hence 
in Laplace space only differ by their frequency depen- 
dent part. They transform into one another by break- 
ing time- ordering between the different parts of the con- 
tour, i.e., freely shifting around vertices without break- 
ing time-ordering on each part of the contour. It follows 
from the diagram rules that they all come with the same 
TMEs, electrode distribution functions, and overall sign, 
cf. Sec. Mil By considering only the time-dependent part 
and its Laplace transform, we derive in Sec. IV Al new di- 
agram rules for evaluating an entire subgroup at once, 
arriving at an expression as simple as that for a single 
diagram. This halves the number of diagrams one needs 
to evaluate - and actually the number can even be halved 
once more. How to achieve this is explained in Sec. IV Bl 
One can exploit that each diagram within a class is re- 
lated to its horizontal neighbor in Fig. [5] by moving the 
latest vertex up or down. Together with a diagram-based 
(rather than rate-based) looping this results in a speedup 
by a factor of 20 in actual numerical calculations. Finally, 
in Sec. IV Cl we explain in what way the classes contribute 
to effective rates for different physical processes and il- 
lustrate the importance of this in the various transport 
regimes defined by the applied bias and gate voltage. 



A. Evaluating subgroups of diagrams 

We start the efficient evaluation of the sum of diagrams 
of a subgroup, G.(x)(t), by selecting a representative di- 
agram and labeling the times r k in all diagrams in the 
subgroup based on this diagram. Here we take the top- 



most diagram in Fig. [5] in each subgroup G.(x)(t), where 
all the vertices on the upper part of the contour are posi- 
tioned at the latest possible times (as far as the subgroup 
allows for this). This choice is advantageous for the fur- 
ther developments in Sec. IV Bl We read off the energy 
differences 8 k from this representative diagram only, and 
Laplace transform the time-dependent part of this indi- 
vidual fourth order diagram [cf. Eq. (|4(il) ]: 
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The other diagrams within the subgroup G.(x)(t) are re- 
lated by breaking the time-ordering between vertices on 
different parts of the contour, but keeping the position 
of the latest vertex fixed at time t. For our choice of the 
representative diagram, this is equivalent to letting the 
vertex with time r 2 move freely. Summing over the three 
diagrams then exactly corresponds to freely integrating 
over r 2 , as is shown in App.[C] Thus the zero-frequency 
Laplace transform of the sum of all diagrams within a 
subgroup G.(x)(t) is given by 

G.(x)(t) ~ / dr 2 / dridr J[ e i(^+i-r k )5 k+1 

J —oo J i n 

*>Ti>T>£ 
p OO 2 p rxj 

= / df 3 e-^' s ^ J[ / df k e-* Skfk 
Jo k=1 Jo 

03 - 02 + «0 d k 

This result is just as simple as for a single diagram; the 
only difference between Eq. (|53a[) and Eq. (|52")l is the en- 
ergy appearing in the leftmost denominator. This is a 
central result of the paper: we directly obtain the con- 
tribution of a whole subgroup by modifying the diagram 
rule for the zero-frequency propagator. One has to eval- 
uate only the representative diagram and assign to the 
latest segment the propagator (S3 — o^ + zO) -1 (instead of 
the usual 5^ 1 ) . This simplification only works under two 
conditions: (i) we are in the zero- frequency limit z — > iO 
and (ii) all secular states are degenerate in energy: either 
E a -E a > > hT (non-secular) or E a -E a > <C hT (secular = 
degenerate) holds, i.e., C ss can be set to zero. In App. [C] 
we show in detail how these conditions enter, in particu- 
lar the proper handling of imaginary convergence factors 
iO, and we discuss a worked out example for subgroup 
C.(x)(t). 

A point which still requires separate care is the secular 
cases of the reducible diagrams in classes B and C. When 
integrating freely over t 2 we sum over all diagrams, in- 
cluding the reducible ones. As discussed in Sec. IIVI this 
should only be done when the intermediate states on the 
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free propagator part are non-secular, i.e., when 82 i= iO 
in B.(l/2)(t) and when S 1 + S 3 ^ S 2 + iO in C.(l/2)(t). 
For intermediate dot states for which this condition is 
not satisfied we must sum up only the irreducible contri- 
butions. However, similarly to the non-secular case, this 
can be effected by a non time-ordered integration over 
r%. For B.(l)(t) and B.(2)(t) two irreducible diagrams 
remain to be summed for the secular case 8 2 = *0: 
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FIG. 6: Example of the construction of gain and loss partners 
in second order. From the diagram rules it follows that the 
contribution of both diagrams are identical except for their 
sign. 
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(53b) 



The modified diagram rule in this case requires 8% for 
the center propagator (instead of the usual S^ 1 )- Note 
that the energies 8k, k = 1, 2, 3, are those of the reducible 
representative diagram, which is actually excluded from 
the sum. For groups C.(l)(t) and C.(2)(t) only the ir- 
reducible representative diagram remains in the secular 
case 82 = 81 + <5 3 — iO, and the standard rule gives 



C.0r)(t) 



8 3 S 2 8i 



(53c) 



The results (|53a[) - (|53c[) can alternatively be obtained by 
directly summing the Laplace transformed propagators 
of the individual diagrams. This is shown in App. [C] us- 
ing general relations between the energy denominators of 
diagrams within a subgroup in the zero- frequency limit. 
This could be of use in diagrammatic calculations of 
quantities other than the density matrix and the current, 
e.g. current noised and time-dependent observables^. 
Finally, we note that for analytic calculations one can 
further sum up the four contributions from each group 
G.(l) to a single expression as well, which is, however, 
of no further advantage for the numerical implementa- 
tion envisaged here, since one exploits the relations (I55al) , 
()55b|) as explained in subsection IV Bl 



The central result Eq. (|53al) can be generalized to any 
order of perturbation theory n, resulting in a relative 
computational gain which grows with n (see App. [C|l . 
By the same three step procedure as outlined for the 
fourth order, the diagrams can be combined into sub- 
groups with x vertices on the upper part of the contour 
(x = 0, . . . , n/2) and y — n — x on the lower one. All 
diagrams in the subgroup are generated by moving ver- 
tices around on the upper and lower part of the contour, 
while keeping the contractions and the vertex at the lat- 
est time t = T n -\ fixed. We sum over all diagrams in the 
subgroup by breaking the relative time ordering of the 



vertices on the different parts of the contour: 



dT„_ 2 • • • d7 



dr y -i • • • dr 



t>T„- 



2"->T y > — OO t>T y 

n-1 



n 

k=y+l 



df fc e" 



-1'">T> — OO 



-(S k -Sy)f k 



n< 

fc=0 
V 

n 



,j-(T k + 1 -T k )8 k + 1 



<\?j(: ' ' " 



n 

k=y+l 



iO 



y 1 



(54) 



The subgroup can thus be summed by using the follow- 
ing modified diagram rule. Determine the propagators 
for the representative diagram as usual. Subtract the 
energy difference 8 V of segment y from all later ones, 
8k — > 8k — S y , k > y. Here 8 y belongs to the segment 
separating the part of the representative diagram with 
vertices only on the upper and lower part of the con- 
tour respectively (ignoring the fixed latest vertex). The 
systematic exclusion of reducible diagrams with secular 
intermediate states (cf. discussion of class B and C in 
fourth order above) can be done most easily in Laplace 
space by extending the method presented in App. [C] 



B. Gain- loss relations between diagram groups and 
diagram-group based looping 

We now shortly address another relation between dia- 
grams, which can be exploited to increase the efficiency 
of their evaluation by another factor two. Each diagram 
in Fig. [5] has a horizontal neighbor, which is identical ex- 
cept for having the last vertex on the opposite part of 
the contour. We will refer to such a pair of diagrams as 
gain-loss-partners, for reasons which become clear in the 
following. Considering a fixed set of intermediate states 
and assuming the final states to be secular (i.e., either 
the same or energetically degenerate) , it is easily verified 
from the diagram rules (see Sec. IB 2\ that moving the 
latest vertex to the opposite part of the contour gives 
the same analytic expression, but with the opposite sign. 
This is illustrated in Fig. [6] This property of pairs of di- 
agrams implies the sum rule (|20[) for the kernel^, (which 
guarantees probability conservation of the density ma- 
trix), but is not equivalent to it. In second order for 
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FIG. 7: Example of the gain-loss-chain among class C dia- 
grams in fourth order. 
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diagonal diagrams as in Fig. |6] the property has a simple 
intuitive interpretation: a tunneling event which changes 
the dot state from a state a to a state b increases the rate 
by which the occupation probability of the final state b 
changes, while it decreases the rate of change of occupa- 
tion probability of the initial state a. The rate for gaining 
probability in state b is described by the left diagram in 
Fig. |6l adding to the kernel element Kg£. Its partner, the 
right diagram in Fig. [5J is obtained by moving the lat- 
est vertex to the opposite part of the contour and gives 
the related rate of probability "loss" for state a. Notice 
that it adds to a different kernel element (namely K^) 
exactly the negative of the "gain" contribution. For nu- 
merical calculations this implies that if one simply loops 
over all possible combinations of initial and final states, 
the same quantity is calculated twice as a contribution 
to two different kernel elements. 

This can be avoided: for problems where only diagonal 
kernel elements (K^)^ need to be determined, one has 
merely to calculate G.(0)(s) and G.(2)(t). This enables 
a very efficient evaluation as follows. For each diagram 
class we take the G.(0)(s) diagram and specify an initial 
state a = a', thereby fixing the final state b = b' = a as 
well (Fig. [7J left most diagram). We then determine all 
allowed intermediate states c\ , C2 , C3 on the lower prop- 
agator. For each such possible sequence of states a,Ci, 
the TMEs need to be evaluated only once per class. We 
furthermore have to calculate only two energy dependent 
functions, one for G.(0)(s) and one for G.(2)(t), and then 
use 



G.(l)(s) = -G.(0)(s), 
G.(l)(t) = -G.(2)(t), 



(55a) 
(55b) 



see Fig. [7J The energy dependent contributions times the 
TMEs can now simply be added to the respective kernel 

super matrix elements (JC^')°J, C^iff Jaa 4 1 see Fig. [7] 
Implementation of this scheme, utilizing the grouping, 
gain-loss relations and storage of/looping over non-zero 
TMEs only results in the speedup of about a total factor 
of 20 for the numerical calculations. 



FIG. 8: Stability diagram (dl/dVb as a function of V g and 
Vb) for the Anderson impurity model in a magnetic field, see 
setup in Fig. [T][a) . The dot single-particle energy for spin- 
projection a — 1 1 1 is given by e CT = eo — eaV g ± Ez/2. Here 
we set the Zeeman splitting to Ez = 0.2(7, the level offset 
to eo = 0.25J7. We used symmetric tunnel rates associated 
with the leads, r s = F d = T = 0.0004 k B T/h and set the 
thermal energy to ksT = 0.01 U. We assume equal capaci- 
tances associated with the source and drain tunnel junctions 
and apply the bias symmetrically. The Coulomb blockade 
regime, where the level is singly occupied, is the central tri- 
angular region within the shown gate range. The positions of 
the second order transport resonances are marked by dashed 
gray lines, the fourth order transport resonances (emphasized 
by white dotted lines) are (i) inelastic cotunneling, (ii) pair 
tunneling and (iii) cotunneling assisted sequential tunneling 
(weakly seen continuation of the Zeeman- lifted | \) — > 1 0/2) 
excitation lines inside the Coulomb blockade region). 



C. Contributions of diagram classes to physical 
processes 

We now illustrate the importance of summing all di- 
agrams of a given order in perturbation theory, and the 
relative influence of the groups of diagrams. This is rel- 
evant for understanding the impact of approximations, 
as well as the relation to the rates calculated in the TM 
approach in Sec. IVI1 For this the simplest models of an 
interacting quantum dot, the Anderson impurity model, 
suffices. This model describes a single level which can 
be populated by at most two (interacting) electrons of 
opposite spin: 



H 



where n a — d\d a is the occupation number and U the 
strength of the on-site Coulomb interaction (excess en- 
ergy required for double occupation). The four many- 
body energy eigenstates are |0) (empty level), \a) (singly 
occupied levels) and |2) (doubly occupied level), with en- 
ergies E = 0, E a = e a and E 2 = J2a ( ° + U ■ Under 
the influence of a magnetic field, the spin-degeneracy is 
lifted by the Zeeman splitting Ez = £f — ej,. 



Fig. [8] shows the corresponding stability diagram, i.e., 
conductance dl/dV^ plotted as function of V g and 14, 
resulting from the full calculation including all second 
and fourth order contributions to the transport kernels. 
We focus on gate voltages around the Coulomb block- 
ade region where the dot is singly occupied. In the cho- 
sen gate range the plot is left-right symmetric (particle- 
hole symmetry) and due to the additional symmetry with 
respect to bias inversion (source-drain symmetry), only 
positive bias voltages are shown. The ground-state to 
ground-state transitions, determining the edges of the 
singly occupied Coulomb blockade region, are due to the 
single-electron tunneling (SET) processes t) ~~ > |0) an d 
|t) — > 1 2). SET transitions involving the spin-excited 
state | l) appear as lines which are separated from the 
ground state transition lines |0) — > t) an d |2) — > t) by 
the Zeeman energy Ez- 

If the kernels were only calculated up to second order, 
only the above mentioned transport resonances would be 
seen in the stability diagram. Fourth order processes give 
rise to three additional types of transport resonances, see 
the numeration in Fig. |SJ 

(i) The horizontal step inside the Coulomb blockade 
region corresponds to the onset of inelastic cotunnel- 
ing&Z&: when the bias voltage exceeds the spin-splitting, 
eV[, > Ez, a coherent tunnel process can take place which 
transfers an electron from the source to the drain elec- 
trode, leaving the dot in the excited |4-)-state. This pro- 
cess involves only virtual occupation of the energetically 
forbidden states |2) and |0) and is therefore only alge- 
braically suppressed by the energy of these states. How- 
ever, as it involves two coherent tunnel processes, it is 
proportional to the fourth power of the tunnel Hamilto- 
nian Ht- Since the charge on the dot is the same before 
and after the cotunncling process, the resonance position 
is independent of V g . 

(ii) There are additional steps in the differential con- 
ductance inside the SET regime, which have the same 
gate dependence as the SET resonances (color change 
along lines ending at the upper figure corners). These 
pair tunneling resonances^ correspond to direct transi- 
tions between the states |0) and |2), involving coherent 
tunneling of an electron pair onto / out off the dot. Note 
that this becomes energetically allowed at a lower bias 
voltage than the sequential addition / removal of two 
electrons (|0) -> \a) —¥ |2) / |2) -> \a) -> |0)). 

(iii) Finally, there are also gate-dependent peaks inside 
the Coulomb blockade regime, the so-called cotunncling 
assisted sequential tunneling (CAST or CO-SET) reso- 
nances^^. These correspond to SET transitions, where 
the initial state is the excited | J,)-state. In second order 
only the ground state is populated inside the Coulomb 
blockade regime, such that this transition cannot take 
place. In fourth order, however, the excited state can be 
populated due to a preceding inelastic cotunneling pro- 
cess, and therefore resonance shows up above the inelas- 
tic cotunncling threshold inside the Coulomb blockade 
regime. 
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FIG. 9: Examples of diagrams with initial states a = a 1 — |t). 
We determine the charge numbers of the allowed intermediate 
and the final states, using that at a vertex where a contraction 
line ends / starts the charge number changes by ±1. Similarly, 
the spin projection changes by ±<r/2 at this vertex, where 
a/2 is the electrode spin-index of the contraction. Using this 
restriction, we find that the group A. (2) can only contribute to 
inelastic cotunneling (b — b' — \D), while C.(2) also allows the 
elastic process (6 = b' = |t)). Groups B.(l) and C.(l) yield 
broadening and level renormalization of sequential tunneling 
processes, like |f) — > |0): B.(l) accounts for the possibility of a 
charge fluctuation in the initial state (|t) — > |0) — > |t)), C.(l) 
for a charge fluctuation in the final state (|0) — > |t) — ► |0)). 

In addition to these resonance effects, fourth order 
terms both broaden and shift the SET resonances and 
give rise to a finite conductance background due to elas- 
tic cotunneling (same as inelastic cotunneling explained 
above, but with initial and final states identical or of the 
same energy). 

We now analyze to which extent the diagram groups 
in Fig. [5] contribute to specific rates for the physical pro- 
cesses mentioned above. In Fig. [§] we illustrate for the 
initial state | f) how the selection rules for charge- and 
spin-projection at each vertex restrict the intermediate 
and the final states. In general, the allowed charge num- 
ber Nb of a final state b for a given initial state a with 
charge N a is readily found for each diagram group G.(x) 
by assigning specific directions to the contraction lines 
and using charge selection rules only. Indicating these 

charge numbers in the kernel by {K^ )jV(,,jv we obtain 
symbolically 

(K { J?) Na , Na ~A.(0) + B.(0) + C.(0) + 

A.(2)+B.(2) + C.(2), (56a) 

~A.(1)+B.(1) + C.(1), (56b) 

N a ±2.N a ~A.(2)+C.(2). (56c) 

From the restricted change in the charge numbers it is 
clear that (|56ap describes cotunncling, (|56b[) corrections 
to SET and (|56c|) pair-tunneling. We note a subtlety for 
(|56a[) : if in addition to the energies, the initial and final 
states are also equal (a = b), the rate must comply with 
the sum-rule (2"0|): 
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FIG. 10: Absolute errors occurring when neglecting contributions of a certain diagram class. Because for class B and C, the 
corrections along the resonance lines (marked by gray dotted lines) exceed the ones in between by orders of magnitude, the color 
scales are chosen logarithmic both in the positive and in the negative regime, i.e. 10 log|(d//dVb) — (d//dVb)approx|. Hereby, 
white color indicates that the contribution from the specific diagram class lies below the threshold for the logarithmic red / 
blue color scale applied for positive / negative error. 



This means that elastic cotunncling rates cannot be 
separated from the "loss" contributions which enforce 
the sum-rule. 

To gain more insight into the physics incorporated in 
the different diagram groups, we now selectively leave 
out contributions and calculate the resulting error in 
the current. By pairwisc neglecting horizontal neigh- 
bors in Fig. [5] i.e., gain-loss partners (cf., Sec. IV B[) . 
the sum-rule (|20|) is conserved. However, the diagram- 
matic grouping reveals that consistency is only guaran- 
teed by neglecting entire classes of contributions. This is 
illustrated by considering the contributions to the pair- 
tunneling rate (|56c|) . Assume that one considers neglect- 
ing the A.(2)(t) contribution in (|56c[) . Then, to preserve 
the sum-rule, we drop in (|56b>|) the contribution of the 
A.(l)(t) subgroup. However, this occurs only in the 
combination A.(l) =A.(l)(t)+A.(l)(s) which contains 
physically necessary^ partial cancellations between the 
two terms: there are pair-tunneling contributions which 
should not influence the SET rate (|56b[) . We therefore 
drop together with A.(l)(t) also the A.(l)(s) term, and 
due to the sum-rule correspondingly in (|56a[) the A.(0)(s) 
term. To keep consistency we thus exclude all diagrams 
G.(2)(t), G.(l)(t), G.(l)(s) andG.(0)(s) of a certain class 
G. 

In Fig. [TUlwe illustrate the impact of the above for the 
Anderson model. Going from left to right we plot the 
absolute error in the differential conductance dl/dVb re- 
sulting from the neglect of cither diagram class A, B or 
C as a whole. Here blue (red) color indicates that in- 
clusion of the specific diagram class reduces (enhances) 
the differential conductance. Lines along which the color 
changes from red to blue indicates an incorrect position of 
the resonance. The occurrence of extended regions with 
uniform color (constant differential error) additionally in- 
dicate that the current is wrong at all voltages above this 



region. 

The left panel in Fig. [TOl reveals that, for the Ander- 
son model, class A does not have any influence below 
the inelastic cotunncling threshold. The reason is that 
by their structure the irreducible diagrams of the group 
A.(2)(t) necessarily involve a spin- flip as shown in Fig. [9]) 
and therefore cannot contribute to the elastic cotunnel- 
ing process |f) — > |f). Class A gives no major correction 
along the SET resonance lines as compared to classes B 
and C below. However, inside in the SET regime the 
increase of the conductance above the pair tunneling res- 
onance shows a large deviation when the class A con- 
tributions to pair tunneling are neglected. Finally, the 
error made in the inelastic cotunncling also affects the 
relaxation of the spin excited state in CAST processes, 
as evidenced by the CAST resonance lines showing up in 
the error plot. 

The situation is completely different when neglecting 
class B diagrams, as shown in the center panel of Fig. [TO] 
Deep inside the Coulomb diamond, they do not give any 
contribution, which is related to their topological struc- 
tured Instead, class B yields significant contributions to 
the SET current level, as exhibited by the uniform posi- 
tive (red) background in the SET regime. Moreover, the 
resonance positions are affected as well: along each en- 
tire resonance line, pronounced "shadows" occur. The 
negative (blue) correction below the resonance and the 
positive correction (red) above signal a shift of the on- 
set of the current, i.e., a level renormalizationd. More 
precisely, class B diagrams can be related to level renor- 
malization of the initial state in a SET process. In SET 
processes, the dot goes from state a to state b by addition 
of an electron, represented by the left diagram in Fig. [5] 
Diagrams from group B.(l) have the same structure, ex- 
cept for an intermediate charge fluctuation ("bubble") 
of the initial state a, see e.g. B.(l)(t) in Fig. [HI with 
a = a' = | t)j b = b' = |0). The lowering of the energy 
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of the initial state shifts the resonance positions of pro- 
cesses where electron tunnels out (in) to lower (higher) 
gate voltages, in agreement with the result in the center 
panel of Fig. [TO] See Refill for a more detailed discus- 
sion. The pair tunneling is not affected since there are 
no B diagrams with two lines connecting upper and lower 
parts of the contour, as required for two-electron transfer. 

In a similar way, class C diagrams account for final 
state level renormalization: the color "shadows" along 
the resonance lines in the right panel of Fig. [10] are prac- 
tically inverted with respect to the result for class B in 
the center panel, indicating an opposite level shift. In- 
deed, diagrams from group C.(l) contribute to tunneling 
from state a to b with an intermediate charge fluctua- 
tion in the final state b, see e.g. C.(l)(s) in Fig. |9] with 
a = a' = | t)> b = b' = |0). Additionally, class C con- 
tributes to inelastic cotunneling and completely mediates 
the elastic cotunneling process | f) — > \ f), see C.(2)(t) 
in Fig. [5] The error at the onset of pair tunneling, to 
which class C contributes, seems less pronounced than 
for class A. In fact the contributions are cqual^i but this 
is masked by the uniform positive correction to SET pro- 
cesses which is larger for class C. 



VI. RELATION TO THE T-MATRIX 
APPROACH 

Since the first studies of higher-order transport pro- 
cesses (see^i), master equations have been used with 
transition rates calculated from a generalized form of 
Fermi's golden rule 3 -! in many works studying transport 
beyond leading order in the tunnel couplin g 22 ' 79 ~ ^2, This 
scattering or TM approach seems to be very similar to 
the GME approaches discussed so far, and in this section, 
we clarify the connection. For this the grouping of the 
kernel contributions discussed in the previous sections is 
of crucial importance. In particular, it reveals the pre- 
cise origin of the divergences occurring in the TM rates 
and allows us to derive the correct regularization of these 
divergences, which differs from the ad-hoc regularization 
employed throughout the literature. 



A. Stationary state equation 

The idea of the TM approach is to apply many-body 
scattering theory in the form of a generalized version of 
Fermi's golden rule to describe stationary quantum trans- 
portal. One calculates the time evolution of the occupa- 
tion probabilities of a state |a) of the total system from 
the transition amplitudes 



b | 5(t) 



(58) 



where T denotes the time-ordering operator. In Eq. (|58[) 
it is assumed that at time t , when the interaction was 
switched on, the total system was in a direct product 



state | a) = |a)|fc) of lead (|&;)) and quantum dot (|a)) 
states. The leads are assumed to be individually at equi- 
librium at that time. This exactly corresponds to the as- 
sumptions of the GME approach discussed in Sec. Ill Al 
As a consequence of the interaction, the state \a) evolves 
into \a(t)) for t > to, which is no longer of product form 
and has an overlap with states \b) ^ \a). The correspond- 
ing transition rate is calculated from the amplitudes via 



_ d 
=Re 



b\a(t) 



at 



b\a(t)))(a(t)\b 



(59) 



Usually, in the TM approach only occupation probabil- 
ities are taken into account, corresponding to diagonal 
elements of the RDM. In general, this can be insuf- 
ficient, as coherences between secular states (degener- 
ate on the scale of T) play an important role for var- 
ious models -e.g. in the case of (pseudo) spin polar- 
izatio n 39 ' 45 . Here we compare the effective rates de- 
termining the occupation probabilities in the TM and 
GME approach. Averaging the transition rate (|59")) . with 
a) = \a)\k) — > \b) = \b)\k'), over the initial (\k)) and 
final states (\k')) of the electrodes with the initial grand- 
canonical probabilities [cf. Eq. (fTTT) ] wc obtain transition 
rates 



Ta^fc = ^2 r r(ak)->(bk>)(k\PR\k) 



(60) 



kk' 



for the time-evolution equation of the RDM, 



Pbb(t) = 2J \Pa->b(t, U)) Paa{to) - T b ^ a (t, t ) Pbbih)] , 
a^b 

(61) 

see^l for a derivation. The rates describe the probability 
for a transition to the state \b) at time t, given that the 
system was prepared in state \a) at time t^. Clearly, in 
the long time limit t — to — > oo this is not an equation for 
the stationary state. Still, the key step in the formulation 
of the TM approach is that one replaces p aa (to) on the 
right-hand-side of Eq. (|61[) by the stationary occupancies 
Paa-, and sets Pbbif) — 0. Then the resulting equation is 
solved for the occupancies, with the rates evaluated in 
the stationary limit to — > — oo. In contrast, in the kinetic 
equation p3[) , the density matrix elements on the right- 
hand-side of the equation are not taken at the initial time 
to, but at times t > to where the system has already 
approached the steady state. We now first show that, as 
a direct result of the above ad-hoc replacement, the TM 
rates calculated to beyond second order in the tunneling 
Ht are divergent in the zero frequency limit, i.e., z — > 0. 



B. Divergence of the stationary TM kernel and its 
proper regularization 

There are two ways to calculate the rates, start- 
ing from the expansion of the time-ordered exponential 
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in Eq. (l58l) : 

[b\a(t) 



1-1 



\-J J drH^(T)j dnff^Ti)---- a). (62) 

To make contact with the literature, we first follow the 
standard route by first performing the time-integrations, 
and obtain from Eq. (|62l) in the stationary limit to — ^ 

— CO 



6 1 5(t) 



= lim 



Er - E a + IT] 



b\T\h 



such that 7 5 _ >b ~ becomes independent of t as expected for 
the stationary state. This gives the well-known general- 
ization of the Golden-Rule rate 3 -^ 



7 .^ = 2v5(E l -E~ a ) 



b\T\ 



(63) 



where the transition amplitude involves the T-matrix T, 
instead of the interaction Ht- The T-matrix is defined 
by a Dyson-like equation 



T = Hn 



H 7 



1 



En - H + ir\ 



-T. 



which can be truncated at the desired order. For compar- 
ison with the GME approaches, it is more convenient to 
alternatively calculate 7 5 _ > { ) to fourth order, postponing 
the time integrations. Setting t% — t in Eq. (|58|) 



3 / -t \ 4 j 3 



n/ ^-i^^-x; 

i=l k^" 7 *" 



x (6 



4-7 3 



1+ £U) II/ ^.iBktk-i) 



Since lead and quantum dot states at the initial time are 
not correlated we can now first trace out the leads. One 
can express the TM transition rates between states of the 
dot a, b as the elements 

r a ^b(Mo) = (&|K™(t-to)|a)H|&) ( 64 ) 
of the superoperator 

JC I TM (t-r) = -Tr R C I T (t)C I T (r)p R 

+ Tr R JdT 2 dT 1 £i(t)£!rto)j%(T 1 )jC%(j)p R . (65) 

t>r 2 >ri>r 

Explicitly, this follows by applying Eq. (|C8|) backwards. 
Alternatively, this equation results straightforwardly 



when formulating scattering theorysi in a Liouvillc or 
"tetradic" formalism^. 

We can now make the connection to the kernel ap- 
pearing in the GME, which we discussed in Sec. Ill Bl . 
and compare the effective rate matrices for the occupan- 
cies. Comparing with the BR superoperator in Eq. pip , 
one finds that its second order part matches the one of 
Eq. (f6"5")) exactly. Going to the Schrodingcr picture and 
taking the Laplace transform [Eq. (|16[) ]. and considering 
super-matrix elements between diagonal states, 



4to 
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Mi 



(66) 



Thus, to the lowest order of perturbation theory, the TM 
approach produces exactly the stationary GME equation 
with a kernel that is well behaved in the stationary limit 
z — > iO. 

However, the fourth order part in (|65p is lacking the 
second term in (|3T|) which subtracts all the reducible 
parts among the fourth order contributions. The physical 
origin of the appearance of the reducible correction term 
was traced clearly in the BR derivation of the GME ker- 
nel in Sec. Ill Bl (and the NZ approach in App. [SJ. There 
the subtraction of reducible contributions, which is miss- 
ing in (|64p . emerged by consistently eliminating p(to) in 
the second order term in favor of p(r). (We note that 
in the RT approach this identification is harder to make, 
since one always deals with correctly regularized expres- 
sions from the start.) This term thus accounts for the 
fact that at times t > to the total system density matrix 
does not factorize anymore, an effect which is however 
only important when going beyond the lowest order. In- 
deed, one directly arrives at the TM approach by ignoring 
this fact in the derivation of the BR approach. Thus, in 
the TM approach one effectively (but tacitly) makes the 
assumption that the dot and reservoir states are statis- 
tically independent after the interaction is switched on. 
We now show that as a result of this assumption the 
rates in the TM approach diverge in the stationary limit. 
Writing out the relation between the fourth order parts 
of TM and BR kernel, it reads in contrast to Eq. 



hi 



K$(z)+K$\z) + K$\z) 



iii, 



(67) 



arises from the second reducible 



fourth order term in (|3ip . which we decompose into two 
parts. The first part contains only non-degenerate (non- 
secular) intermediate states, 



( ^nn) 



K£}{z). 



(68) 



This is precisely the correction term to Kss in the effec- 
tive master equation for the secular part of the density 
matrix (which here reduces to the diagonal part) arising 
from the non-diagonal coherences, cf. Eq. (|^|) . As ex- 
plained in Sec. II VI it must be included to obtain a system- 
atic expansion of the effective transition rates between 
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probabilities in powers of T. Both kernels, Kii\z) and 
K^\z), are well-behaved for z —> iO (see Sec.|VJ). The 
remaining term contains only intermediate states which 
are strictly degenerate in energy (strictly secular) 



K, 



Pi*) 



K 



(2) 



z 



(69) 



and as a result diverges as z 1 since Kss(z) is well- 
behaved for z -> iO, cf. Eq. ([66]). Therefore the TM rate 

K-tm( z ) diverges as z _1 as well. Rewriting Eq. ([67]) we 
can express the effective GME kernel P9"]l determining 
the stationary occupation probabilities: 



(70a) 
(70b) 



This equation summarizes the central relation of GME 
approaches to an automatically regularized TM expres- 
sion. The effective fourth order kernel for the proba- 
bilities is thus obtained by either in Eq. (|70ap adding 
the non-secular reducible correction to the GME kernel 
(both finite) or in Eq. (|70b[) subtracting from the TM 
kernel (|55t the secular reducible correction and canceling 
the z~ l divergences. We have thereby precisely identi- 
fied the correct regularizing term (|69[) which should be 
used if one would like to keep on using a fourth order 
TM, expressed in terms of known second order TM rate 
expressions [cf. Eq. ([6"6"]) ] and the frequency z. We em- 
phasize furthermore that the form (|70a.|) makes explicit 
that the probabilities contain corrections from the non- 
secular coherences whereas (|70b[) and, in fact, the TM 
approach itself, make no reference to non-diagonal den- 
sity matrix elements. 

Eq. (|70b[) shows explicitly that the relevant kernel for 
transport and scattering problem are related. General 
relations between irreducible kernels and Liouville TM 
expressions were given first by Fano^. However, to our 
knowledge, the relation between the TM rates and the 
non-secular corrections to the effective fourth order ker- 
nel, has not been addressed before. 



C. Regularization error for Anderson model 

A key point of the TM approach as formulated in the 
literature is that instead of Eq. (|70b[) one uses an ad- 
hoc regularization to cure the divergence which was in- 
advertently introduced by the ad-hoc formulation of the 
stationary state equation. The resulting TM rates for 
occupations obtained in the literature show a striking 
similarity to the GME expressions for the effective rates 
for the occupancies, including the non-secular contribu- 
tions from the coherences. The reason for the similarity 
is also made explicit in Eq. (|70b[) . However, we now ex- 
plicitly show that the ad-hoc regularization significantly 
differs from the correct regularizing term Eq. (|69l) sub- 
tracted in Eq. (|70b[) . To illustrate its importance ana- 
lytically, we consider the simple Anderson model already 



studied in Sec. IV CI In applications of the TM approach, 
typically not all fourth order contributions are included. 
Very commonly, corrections to SET, i.e., diagram groups 
A.(l), C.(l) and the complete class B, are dropped a 
priori. As this breaks the gain- loss-chain (see Sec. IV Cf . 
the sum-rule is then enforced by hand, which makes the 
groups A.(0) and C.(0) redundant as well. 

As we do not wish to study here the errors arising 
from such additional approximations (see Refi22 for such 
a comparison), we focus on a cotunneling contribution 
which would be included in any type of TM based cal- 
culation, namely the one from diagram group C.(2).(£) 
(contributions from A.(2).(i) drop out for infinite charg- 
ing energy) to the elastic process \a) — > |0) — > \a). We 
compare the energy dependent parts of its analytical con- 
tribution to the kernel element K%% for the GME, respec- 
tively (.Ktm)^ f° r the TM. In Sec. [Vj we determined the 
energy dependent part of the contribution of subgroup 
C.(2)(t) to be [Eqs. (|55a|) and ([55c])]: 



C J Pi *a (*3 - + *0)] 1 S 2 / 6i + 6 3 - iO, 

' [ ,U ~ I Mafc] -1 S 2 = 6! + S 3 - iO. 

This function clearly distinguishes between the cases of 
secular and non-secular intermediate states, avoiding the 
inclusion of the divergent reducible term. For our ex- 
ample, we determine 81,82,83 diagrammatically using 
Fig. 0] For the kernel element in the Anderson 

model, the states are a = a' = b = b' = |cr), c = d = |0), 
where a 6 {t,4-}: 

So = 

£1 = iO + E a - E a +J, 

82 = iO — uj + u>' , 

S 3 = iO + E a -E -u). 

The corresponding TM expression misses the exclusion 
of the reducible diagrams and is always given by 



C™(2)(t) 



[Si S 2 (S3 - 6 2 

1 



i or 1 



1 



uj - uj 1 - iO |u/ + e q - E a + i0\ 



even when the secular condition 5 2 +i0 = 8i+S 3 is fulfilled 
here. In the final expression for the rate, C™(x)(t) is 
multiplied by two Fermi functions and we integrate over 
their arguments, uj and u>' respectively. Performing the ui 1 
integral first, one is left with a divergent ui integral due to 
the modulus-square factor in C™(2)(t). The standard 
way to regularize the TM rates mentioned above now 
proceeds as follow o 80 i 82 . For any function F(ui) which is 
well-behaved at ui = one expands in the infinitesimal 
iO 



duj- 



F(u) 



= / dw- 



F(0) 



(ioy 



duj 



F{u>) 



F(uj) - F(0) 
uj 2 - (iO) 2 

(71) 



(iO) 2 1 
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where j' denotes a principal part integration. Divergent 
contributions oc 1/ (iO) are claimed to be due to sequential 
tunnel processes which are already included in other rates 
and are ignored. As shown above in general, however, 
the divergent term is due to neglecting the mixing of the 
lead and dot states, which is not an effect of sequential 
tunneling (see Eq. (|66|) h Moreover, the regularization 
procedure docs in fact not reproduce the regularization 
which is automatically included in the GME approach. 



To illustrate this quantitatively and gain insight into 
which voltage regimes this matters, we have plotted in 
Fig. [TT] (top panel) the differential conductance com- 
puted with the TM approach, using the same model 
and parameters as in Fig. [8] To demonstrate that ex- 
isting discrepancies are not simply healed by including 
the typically neglected corrections to SET, we have em- 
ployed here the "best-possible" TM kernel, taking into 
account all contributions arising from Eq. (|65[) . just reg- 
ularizing the occurring divergences according to Eq. ([71]) . 
For diagram class B, such regularization leads indeed to 
a complete omission when 62 = iO (secular intermediate 
free propagating states), which is the case for the An- 
derson model in our example. Class A, containing only 
irreducible diagrams, requires no regularization and is 
taken into account fully. Both the GME kernel as well 
as this "best-possible" TM kernel satisfy the probability- 
conservation sum-rule. The associated current kernels 
are constructed as discussed in Sec. Ill A[ 



In the lower panel we show the relative deviation be- 
tween the TM and GME results, blue (red) color indi- 
cating that the exact (i.e. GME) differential conduc- 
tance falls below (exceeds) the one obtained from TM ap- 
proach. Clearly, the agreement is good only deep inside 
the Coulomb blockade region. However, all resonance 
lines are dressed by a pronounced red (blue) shadow from 
above (below) , indicating that the corrections to sequen- 
tial tunneling (level renormalization, broadening) are not 
correctly taken into account in the TM. For this simple 
setup, the maximum of deviation encountered amounts 
to up 30% overestimation and 5% underestimation of the 
correct result. Given the parameters of the model, we 
have actually given with Fig. [TT] the best result that can 
be possibly obtained within the TM approach with the 
Eq. (fTTj) . This includes the effects of clastic and inelas- 
tic cotunncling, pair tunneling and single-electron level 
renormalization and broadening effects. Often other sim- 
plifications are made, in addition to the above procedure 
employed, leading to further deviations. For instance, 
quite commonly, only the fourth order cotunneling rates 
from the diagram groups A. (2) and C.(2) are taken into 
account. The figure shows that the TM approach basi- 
cally only works in the "deep Coulomb blockade regime" 
where SET and CO-SET processes are suppressed. 
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FIG. 11: Top: The stability diagram for a single-level quan- 
tum dot as in Fig. [8] calculated using the TM approach. Mid- 
dle: Absolute deviation between the non-linear conductance 
calculated with the GME and TM approach, in logarithmic 
scale, i.e. 10 log|(d/GME/dVb) — (d7rM/dVb)|, with color coding 
just as in Fig. 1101 Bottom: Relative deviation between the 
non-linear conductance calculated with the GME and TM 
approach, [(d/GME/dVb) - (d7 T M)/dV b )] /(d7 G ME/dF b ). Al- 
though the upper panel seems similar to the GME result in 
Fig. [5] indeed large relative errors exist in the vicinity of res- 
onances and throughout the regime where transport is not 
suppressed by Coulomb blockade, in particular when electron 
pair-tunneling in the single-electron tunneling regime^ is en- 
ergetically allowed. 
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VII. SUMMARY 

In summary, we have studied the systematic calcula- 
tion of the fourth order kernel within the generalized 
master equation approach for transport through quan- 
tum dots. At present this is the only approach which 
can efficiently deal with strongly interacting systems with 
complex excitation spectra in the non-linear transport 
regime, while accessing the regime of moderate tunnel 
coupling by a perturbative treatment of higher order tun- 
nel processes. The precise understanding of the calcu- 
lation of the kernel determining the transport rates is 
therefore of great practical importance. On the one hand, 
simplifications which speed up numerical calculations are 
crucial to allow more complex physics to be addressed. 
On the other hand, the comparison of calculations with 
kernels evaluated using different methods is an urgent 
issue. We summarize our main achievements: 

We have first shown the equivalence between the 
real-time diagrammatic approach (RT) and the Bloch- 
Rcdficld quantum master equation (BR) (and the 
Nakajima-Zwanzig technique). This was done both for- 
mally and explicitly, by mapping irreducible RT diagrams 
onto operator expressions in the BR. In particular, we 
showed that in the BR, the well-behaved kernel from the 
diagrammatic approach is obtained as a sum of two parts: 
One part contains all (irreducible and reducible) terms 
of a certain order. The other part exactly cancels all the 
reducible terms, which diverge in the stationary (zero- 
frequency) limit. 

Next we addressed the calculation of the density ma- 
trix using this kernel. The commonly used secular ap- 
proximation for the reduced density matrix was previ- 
ously shown to break down when going to next-to-leading 
order in the tunneling^. Despite this, an effective kernel 
determining the secular part of the reduced density ma- 
trix can be derived: Here we showed that the non-secular 
corrections can effectively be accounted for with a secular 
kernel by including certain reducible diagrams with non- 
secular intermediate states. By adding these to the stan- 
dard irreducible diagrams a new structure is revealed. All 
diagrams can be sorted into groups of diagrams with the 
same distribution of vertices over the Keldysh contour 
and the same contractions. We showed that subgroups 
consisting of three diagrams can be summed analytically, 
thereby avoiding unnecessary integral evaluations. We 
derived new diagram rules for evaluating an entire sub- 
group at once, arriving at an expression as simple as that 
for a single diagram. The physics behind diagrams and 
groups was illuminated by a study of specific contribu- 
tions for a single level with finite Coulomb interaction 
and Zceman splitting (single impurity Anderson model). 

Furthermore, the summation of subgroups allows for a 
general comparison between the GME and TM approach 
("generalized Golden Rule"): We first showed formally 
that the TM "rate-kernel" equals the correct fourth or- 
der kernel plus a divergent term. Thereby we identified 
precisely the correct term needed to regularize the TM 



rates. We emphasize that the GME (BR or RT) kernels 
are well-behaved and finite by construction, the regular- 
ization is incorporated automatically and not put in "by 
hand" . We showed both numerically and analytically 
that the regularizations existing in the literature are in- 
correct and lead to analytical deviations in the individ- 
ual rates and, as a result, to pronounced errors in the 
calculated transport current if one is not deep inside the 
Coulomb blockade regime. We illustrated this numeri- 
cally for the example of the Anderson model in magnetic 
field. 
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Appendix A: Nakajima-Zwanzig projection 
technique 

The standard Nakajima-Zwanzig projection operator 
techniqu o 58 ' allows for a compact and concise deriva- 
tion of an exact expression for the kernel, see e.gi 60 ' 61 and 
the references therein. We briefly review the notation and 
the steps involved in the derivation of the time evolution 
kernel in the time-domain which is performed concisely 
in the interaction picture. In a similar fashion, one can 
derive the current kernel. Using projectors V = puT^R 
and Q = 1 — V we decompose the total density matrix, 
p[ ot = Vp[ ot + Qp( ot , and project the interaction-picture 
Liouville equation (|24|) for the full system: 

Vp{ ot = -iV£\{t)Q P { ot , (AI) 
Q P { ot = -i Q-Cj(i) Qp( ot - i QC{(t)Vp( 0V (A2) 

Here, the crucial property VL\(t)V = was used, which 
is due to the fact that the tunnel Hamiltonian ©, and 
thus also £y ([7]), contains exactly one lead operator: the 
trace must yield zero. Next, the second equation is for- 
mally integrated using Qptot(^o) — and treating the 
term with Vp( ot as a given inhomogencous term: 

Sptot(t) = -i f&T.Te-'^^^QCk^VpUTi)- 

Jta 

(A3) 

Here T is the time-ordering superoperator. Substitution 
into Eq. (|A1|) gives the kinetic equation (|30|) with the 
formally exact kernel: 

/C'(t,r) = — Tr r (T^(i)e-^ d -^)^(r)p R ) , 

(A4) 
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which transforms into the well-known result for the time 
evolution kernel in the Schrodinger picture, Eq. ([T4"|). 
The time evolution kernel (|A4|) contains the non-trivial 
evolution operator which can be expanded in the pertur- 
bation QC^t^Q, 

e -*a/*dn£4(n)Q = 1 _ ; | ' dri Q£{( Tl )Q 

- Jd^dn Q£ i t (t 2 )Q£ i t {t 1 )Q + ... . (A5) 

t>T2>Tl >r 

Applying this we obtain the kernel to fourth order, which 
can be compared to Eq. pil) to confirm the equivalence 
to the BR approach. Inserting Q = 1 — V everywhere 
and using that V gives non-zero only when acting on an 
even number of Cj by Wick's theorem we obtain 

/C J (i,r) = -Tr R C{(t)C{(r)p R + 

Tr r J dr 2 dri£((i)4(r 2 )(l-P)4(T 1 )4(r)p R . (A6) 

t>T2>T!>T 

The first fourth order term, involving the unit operator 1 
in the middle, gives rise to all possible contractions from 
which the reducible ones are subtracted by the second 
term with V in the middle. Alternatively, one may also 
first integrate Eq. (|A1[) for V p\ ot with initial condition 
Vp I tot (t ) = p I tot (t ): 

VpUt) = Ptot(to) - i f dT 2 7^(T 2 )Qp t 7 ot (r 2 ). (A7) 
Jt„ 



and fourth order kernel contributions for the Anderson 
model, starting from the expressions given in Fig. [3J for 
readers not familiar with either the RT or BR technique. 
In App. IB 2l we summarize the general diagram rules, i.e., 
how to set up the time evolution kernel in terms of di- 
agrams and how to directly read off the final expression 
for the contribution from each diagram. 



1. Examples 

We consider the Anderson impurity model with finite 
Coulomb interaction introduced in Sec. IV CI character- 
ized by the four many-body states |0), |t), |2), 
corresponding respectively to zero, one spin-up, one 
spin-down or two spin-paired electrons on the dot. We 
demonstrate the technique by evaluating only the kernel 
element {K)\l = (2| [if|2)(2|] |2), which contains in 
second order the two "loss" rates for processes |2) — > \<r). 
In fourth order it includes besides "loss" rates also the 
elastic cotunneling |2) — > |2). 

For the calculations we will need the following, gener- 
ally valid transformation for kernel elements from time 
space to Laplace space, 



-,-lHt 



-**(t-T')] e *** b'), (Bl) 



Substitution of Eq. (|A3j) into the right hand side of 
Eq. (|A7[) and taking the trace gives a Dyson-type integro- 
differential equation for the reduced density operator 
p 1 = TrR^Ptot- The equivalent equation for the prop- 
agator defined by p 1 (t) = TT 1 (t, io)p J (io) reads 

7T / (Mo) = l+ J dT 2 dn)C I {T2,n)7T I (n,to), (A8) 

t>T2>T 1 >t 

with the kernel given by Eq. (|A4|) . revealing total equiv- 
alence to Eq. (|3T|) as obtained in the RT approach. 



Appendix B: Time- and frequency space calculation 
of the kernels 



In this appendix we present the details of the calcula- 
tion of the kernel Eq. (|14p starting from the interaction- 
picture expansion used in the BR approach in Fig. [3] We 
explicitly obtain a result which was mentioned in Sec. IIIII 
and used as a starting point in Sec. [V] of the main text: 
contributions to the kernel represented by diagrams dif- 
fering only by relative time-ordering of vertices on dif- 
ferent parts of the contour, deviate only in the time- 
dependent function or its Laplace transform (|46p . In 
App. IB II we first discuss an example calculation of second 



where t = t — t. As compared to Eq. (fT6|) . additional 
exponentials arise because of the fact that matrix ele- 
ments with respect to the states of the RDM - which 
transforms according to Eq. (j32|) - have been taken. 



2nd order 

From Fig. [3]wc can infer which diagrams contribute by 
using the charge-selection rule: at each vertex the charge 
changes by one. For the Anderson model, starting from 
1 2} only the intermediate state |cr) is possible, with either 
spin a 4-- The zero-frequency contribution reads 
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This is a sum of two complex conjugate expressions, so 
we only need to evaluate e.g. the last complex expres- 
sion. From Fig. [3] we obtain the time-dependent inter- 
action picture expression, which we transform to the 
Schrodinger picture and Laplace transform with respect 
to the time-interval r' spanned by the diagram and send 
z — > iO 



24 



lim 

:-H0+ 



dr' e izT '(CoC+) 



x(2| e~i HT ' |2)-(2| e-i H ^- T '^D+ \a) (<r| D3 ei Ht |2) . 

Here the expressions left and right of the • correspond to 
the upper and lower contour, respectively. The exponen- 
tials containing H arise from transforming the operator 
expression from the interaction to the Schrodinger pic- 
ture. Next we transform the occurring operators D and 
C to the Schrodinger picture, according to Eqs. ()44a[) and 
(|44b|) . With the use of Eqs. (|4l2"a|) and (|42b]) as well as 
Eqs. (|5aj) and (|5b|) we find 



dw pl S {u) 



where a = -a. Here = f{{uj - Pi)/(k B T)) = 

( e («-w)/(*BT) + i s the Fermi function of lead I 

with temperature T, = 1 — fi{oS), and pia(u>) 



a, i 



as it occurs in Eq. ([5]), is the (possi- 
bly spin-dependent) density of states in lead I. Notice 
that for our simple Anderson model we have no q depen- 
dence of the single particle tunneling amplitudes, and 
thus there is no q dependence of the TMEs here. Fur- 
thermore, the time t cancels out in the exponentials after 
performing the t' integral. This example illustrates the 
form of the time-evolution factor of any 2nd order con- 
tribution: 

1 1 



where Ao is the sum of energies occurring in the argu- 
ment of the exponential which contains the time r'. 



4th order 

For the Anderson model with non-magnetic electrodes, 
selection rules cause the fourth order non-secular correc- 
tions to vanish, i.e., K ns = in Eq. ([50]) . Therefore only 
irreducible contributions remain. These selection rules 
can be used furthermore, in addition to the charge selec- 
tion rule, to determine the allowed intermediate states on 
the diagrams below. For the example kernel element the 
charge number does not change, and therefore all con- 
tributing diagrams have an even number of vertices on 
each contour (i.e., groups G.(0) and G.(2) in Fig. [5]). No- 
tice further that there are no diagrams from group A. (2), 
because these would involve an intermediate charge state 
with three electrons on the dot. 



2« 2 2<e 

" \ 2 



^S^ 2 2^S^ 

a 2 a a 2 a 

/ 2 - " 



2 9»> ,2 2 ^TV 



a a -j- 
a a 



2 2 
+ 

2 2 



+ c.c. 



We calculate the last shown diagram, starting from 
the expression given in Fig. [3] and introduce the time- 
distance of vertex i to the final time t[ := t — Tj, and 
perform the same steps as in the second order example: 



- lim fdr'e- f dn /' dr 2 (C "C+> (C+C£) 

> > 2 *-H0+ J J t _ T , J T1 



x (2 e~i m D 



£>2 e^ H{t - T ' 



S " 4 E F At ' f d < A r 2 / d0J I <^'pro{u)Pi> 
,,, Jo Jo Jo J J 



X/+H /,7(w') T+(2, a)T^(a, 2)T+- (2, &)TjZ{a, 2) 
-iE / dw 1^' |r { +(2, < 7)| 2 |2^(2 )£ r)| 



1 



{-u+Ei-E a )r[ j-(-uj'+E 2 -Ev)T' j.(uj' -E 2 +E^+i0)r' 



-uj + uj' + iO -u> + E 2 - E a + iO u>' - E 2 + E a + iO 



In this example all contracted electrode operators have 
inverted time-order (corresponding to the earliest vertex 
being on the lower contour): as a result the U) contraction 
gets the / ; + function in contrast to the contraction in the 



second order example. We note that the structure of the 
Laplace transformed time-evolution factor appearing in 
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the fourth order contributions has the general form: 
11 1 1 



A A + Ai A + Ai + A 2 8iS 2 S 3 

where Ao/1/2 denote the arguments of the exponentials 
containing the times t'/t^/t^, respectively. 



2. Diagram rules for zero-frequency kernel 

We now give the rules by which one can directly write 
down the diagrammatic representation of the effective 
time evolution kernel and afterwards simply read off from 
each of these single (ungrouped) diagrams the resulting 

analytical contribution to K^\z), Eq. (|49|) . in the zero- 
frequency limit z — > iO. We illustrate how this applies to 
the examples of the previous section. Notice that modi- 
fications of these "traditional" diagram rules in order to 
account for a whole subgroup of diagrams were presented 
in Sec. 



The kernel in diagrammatic representation 

The rules for drawing all diagrams representing the 
effective kernel of even order n = 2,4, .. are as follows: 

• Draw all distributions of n vertices over the two 
contours, vertex k being at time Tfe, k = 0, . . . , n — 
1. Vertices n and 1 are at the boundaries of the 
diagram at times r„_i = t and tq = t, respectively. 

• For each distribution, contract all n vertices in 
n/2 pairs, denoting each contraction by a directed 
line. Each resulting diagram represents a distinct 
contribution to the effective time evolution kernel. 
Note that all irreducible and reducible contractions 
need to be included, where irreducible diagrams are 
those which can nowhere be vertically cut without 
cutting a contraction line. 

• To each contraction j (1 < j < n/2), assign an 
energy uij, as well as lead and spin indices lj, (Xj, 
respectively. 

• On each contour, assign to each segment between 
two vertices a many body state of the quantum dot . 

For the previous examples in second and fourth order we 
thus obtain: 

2< 2 



2nd 



4th 




Translating a diagram 

The rules for translating a diagram into an analytical 
expression can be divided into rules for determining three 
factors. 

1. For each contraction j write a Fermi distribution 
function and spin-dependent density of states 



(B2) 



where + (— ) is chosen if the contraction line agrees 
(disagrees) with the contour direction. For the ver- 
tex at which the contraction starts / ends write a 
many-body tunnel matrix element 

where a and b are the states before and after the 
vertex, respectively, following the direction of the 
contour (not of time). 



2nd 



4th 



/r( W )^( W )T+(2,a)T^( ff ,2) 

/i + H fv{u')piB{u)pi>s{uj') 
xT+(2, a)T^ s {a, 2)T+ S {2, a)T^{o, 2) 



2. Determine the time evolution factor 

n-l 



(B3) 



k=0 



by drawing through each segment of the diagram 
between consecutive times Tk and t^+i (0 < k < 
n — 2) a vertical cut (see Figs. El [7] in the main 
part of the text). Obtain the denominator Sk by 
adding / subtracting the energies of the dot states 
on the contour and the energies of the contractions 
depending on whether they hit the vertical cut from 
the left or from the right. Also add to Sk the fre- 
quency z = iO. 



2nd 



Si + iO = E a + uj - E 2 + iO 



4th 



S3 
S 2 
Si 



■iO 
iO 
iO 



E 2 
E 2 
E a 



uj + E a + iO, 
lj + uj' - E 2 ~ 
-U/-E2 



iO, 



iO. 



(Note that iO is not a convergence factor put in by 
hand, but naturally arises from the Laplace trans- 
form, reflecting the correct analytic behavior of the 
kernel.) 

3. For the diagram as a whole determine the phase 



h 



(B4) 



by counting the number of crossing contraction 
lines n c , and the number of vertices on the lower 
contour n\. 

, _±f _i\0+2 _ _i_ 



2nd 4th 



2G 



Finally one multiplies the three factors, integrates over 
all frequencies uij and sums over all spin values o~j and 
electrodes lj. Notice that also all possibilities for inter- 
mediate quantum dot many body states on the contour 
have to be summed over. The diagram rules, as formu- 
lated above, provide the key insight needed in the main 
text: since diagrams within a group, as defined in Sec. |V1 
are related by moving their vertices around on each part 
of the contour, only the factor arising from the time- 
evolution is different. The ordering of the vertices, the 
direction of the contractions relative to the contour, as 
well as the number of crossing lines and vertices are all 
preserved under this operation. 

The derivation of the rules in the form presented above 
can be found in^. Because of the importance of the time- 
evolution factors we comment on the relation between 5k 
as defined by 



8 3 d 2 5 1 S 



fe+i 

1=0 



(B5) 



and the value obtained by diagram rule 2. This is easily 
seen once one notices the diagrammatic meaning of A; 
(see£2, p. 95): it equals the sum of energies of all lines 
going into the vertex I minus those of the lines going out, 
i.e., the dot energies of the in- and out-going lines of the 
contour and the energy of the contraction which starts / 
ends at vertex I. Summing contributions from all vertices 
from earlier times 77 < Tk , the energies of all contractions 
which have started and ended thus cancel out, leaving the 
difference of energies of contractions running backward 
and forward with respect to time. In the sum, the dot 
energies of subsequent vertices on the same part of the 
contour cancel out, leaving only the difference between 
the upper and lower contour dot energy between vertex 
k and k — 1. Summing all from both contours, one 
gets X)fc=o Afc = This condition ensures that finally no 
exponential containing t is left, reflecting that the kernel, 
see Eq. ()13[) . depends merely on the time difference r' = 
t — t. A more explicit version of this proof can be found 
in Ref3. 



Appendix C: Diagram grouping 

In this appendix we prove the statement made in 
Sec. El freely integrating over the intermediate time r 2 
in the representative diagrams is equivalent to summing 
the three diagrams within the corresponding subgroup, 
provided that the initial states (earliest times) are degen- 
erate. Additionally, we show how to easily calculate the 
partial summation of one or two irreducible diagrams in 
a subgroup in the case where secular diagrams have to 
be excluded. This exclusion is automatically obtained 
in our transport theory and is an important result of 
the paper: it prevents the divergences which plague the 
TM approach from appearing. We present the deriva- 
tion both in time- and frequency-representation which 



FIG. 12: From the representative diagram, all members of the 
triple group can be constructed, exemplified here by subgroup 
C.(2)(t). 



each have their distinct advantages. 

As explained in the main text, we want to express the 
sum of diagrams in a subgroup in terms of the contri- 
bution of a representative diagram (topmost diagram in 
each subgroup in Fig. [5|. The standard diagram rules ex- 
press this contribution as the product of propagators over 
time- intervals of length —fk = Tk — Tk-i, see Eq. (f52"| : 



- j- (5 3 f 3 +S 2 f 2 +Sl Tl ) 



G.(z)(t) ~ / dndf 2 df 3 e 



dfidf 2 df3 g(T 3 ,T2,Ti,To). (CI) 



Here 5k is the sum of energies of the backward minus the 
forward moving contour parts and contraction lines in the 
segment of the diagram between times Tk and Tk-i, see 
Fig. HI including also the Laplace variable z — > iO. The 
simplifications only work for the zero-frequency Laplace 
transform of the kernel which is all that is needed for the 
stationary state. 

We now allow the next-to-last vertex (at time t 2 ) in the 
representative diagrams to move to earlier times (to the 
right in the diagram), see Fig. [T^] Thereby we generate 
the other diagrams in the subgroup. In the first step, 
the vertex at t 2 is permuted with the one at time t\ and 
thereby the energy difference Si of the segment bounded 
by these vertices is changed to a different value: <5 2 — > 5' 2 . 
Permuting the vertex (now at time t\ ) in the second step 
with the vertex at To, the energy difference Si changes as 
well: Si — > S[. From the diagram rules it follows that for 
each such permutation the following relation holds: 



$n + S' n = S' n+ i + S 



(C2) 



i.e. the average of the old and new value (left) equals the 
average of energy differences of the two adjacent segments 
(right). Note that on the left the energy has already been 
modified to S' n+ i by the preceding permutations (with the 
exception of the latest one £3). This is the key relation al- 
lowing the summation of diagrams within each subgroup. 
We now apply this to the three diagrams in Fig. [T2] and 
obtain 



S'o 



Si, 



Si + S'i = S' 2 + (Jo- 
Combining these we obtain the relation 
S-2 + S[ = 83 + So . 



(C3) 
(C4) 

(C5) 
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Here So = E a — E a > denotes the energy difference of the 
initial states, i.e. outside the diagram, and does not in- 
clude the Laplace variable z = iO as the 8k, k = 1,2,3 
do. Using this, the contributions from the generated di- 
agrams can now be expressed in the representative func- 
tion, with the time-arguments of the corresponding ver- 
tices permuted: 

e-^fcft+aSft+fm) = ff (r 3 ,r 1 ,T2 ) ro) ! (C6) 

t t 

e-*( 53f3+5 = f2+ ^ fl ) = 5 (r 3 ,ro,r 2 ,r 1 ). (C7) 

t t 



The arrows indicate how, by permuting time-arguments, 
flCBJ) is obtained from ([CT]) , and (|U7j) from $C6$. Al- 
though this result may seem obvious, one should note 
that the last equation only holds under the condition 
Sq = 0. Thus, the diagram obtained by permuting the 
latest vertex with one of the earlier ones can only be ex- 
pressed in the representative diagram if the initial states 
arc degenerate. Time-ordered integration of the sum of 
all the diagrams in the subgroup is now seen to be equiv- 
alent to decoupled time-integrations on opposite parts of 
the Keldysh contour of the single representative diagram, 
neglecting the remaining diagrams from the group: 



G -(£)(t)~ JdT 2 dridro [g(r 3 , r 2 , n, T )+g(r 3 , n, r 2 , T )+g(r 3 , r , r 2 , n)] = /dridr / dr 2 g(r 3 , r 2 , n, r ) 

T"3>-T2>Ti>T > — OO T 3 >T 1 >T >~ CO T 3 >T 2 > — OO 

(C8) 



Inserting the form of g from Eq. (|C1[) , and changing vari- 
ables to time-intervals on the separate contours, t 2 = 
t 3 — Ti, T\ = ri — To (forward) and r'i = t 3 — r 2 (back- 
ward) the integrals decouple as usual. We obtain the 
result in the main part of the paper: 



G.(x)(t) 



(iO + 8 3 - 82)8281 



non-secular 



(C9) 



Note that the iO has to be supplied "by hand" explicitly 
since it formally cancels in the difference £3 — 5 2 . Below 
we show that this does not alter the value of the integral 
and that the sign automatically follows from the correctly 
regularized terms which arc being summed. 

For diagram class A this completes the derivation. 
However, diagram classes B and C contain reducible dia- 
grams, which diverge if one would allow for secular inter- 
mediate states. In this case, one thus has to perform a 
partial sum of the irreducible diagrams in the subgroup 
only. Note that in our formalism the exclusion of these 
cases is automatically enforced, i.e., we do not exclude 
them "by hand" based on the mere inconvenience of di- 
vergent terms. 

Although the explicit result for classes B and C can 
be obtained in the same way as above, we now show 
that here the frequency space representation has definite 
advantages. 

We first derive Eq. (|C9[) again in frequency-space by 
directly summing the Laplace transforms of the prop- 
agators [left hand side of Eq. (|C1|) and Eqs. (|C6|) and 

1 1 1 



G.(z)(t) 



s 3 S2S 1 836^ Wi 
111 



8382S1 <5 3 ^(5i S'^281 



non-secular 



(CIO) 



We first performed the partial sum over the last two 
terms using relation (|C4|) and assuming 8q = 0. Adding 
the last term and using relation (|C3[) we obtain the full 
sum. Expressing 8[ in the energy differences of the repre- 
sentative diagram 83,82,81 using Eq. (|C5[) we again ob- 
tain Eq. (|C9|) . However, in addition we have treated the 
cases of secular intermediate states as well. For the B.(l) 
and B.(2) subgroups, the representative diagram is itself 
reducible and the secular case arises for <5 2 = iO. If it is 
excluded, we keep only the partial sum in the second line 
of Eq. (ICTOl) : 



B.(x)(t) 



1 



818, 



secular 



(Gil) 



Note that here the iO needs not be written out explic- 
itly. In contrast, for the C.(l) and C.(2) subgroups, the 
representative diagram is the only irreducible one which 
has to be kept: inspecting the reducible diagrams we see 
that the secular case arises for 8' 2 = 83 — 82 + 8\ = iO and 
we can eliminate one parameter. Choosing £3 we obtain 



C.(z)(t) 



1 



£3 (#3 + ^1)^1 



secular 



(C12) 



We note that also here iO needs not be written out explic- 
itly since the 2i0 in £3 + 81 are sufficient to guarantee the 
correct analytic behavior as function of the frequencies. 
We now shortly comment on this point as well on the iO 
explicitly added in Eq. (|C9|) . 

Above we have calculated the time-integrals of the 
time-evolution factors in the diagrams only. The re- 
sulting expressions, multiplied by the statistical factors 
(Fermi- functions) , still need to be integrated over the fre- 
quencies of the contractions which are included in the 
energy differences 8k- These integrals are exactly those 
of the standard perturbation theory and can be found 
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in, e.g., Refs i 27 i 75 . However, in the above (partial) sum- 
mations over subgroup diagrams, we have at several in- 
stances used that we can replace iO by 2i0. This does 
not alter these integrals. The integrand possess a count- 
able number of poles and decays sufficiently fast for the 
residue theorem to apply. Closing the integration con- 
tour in the upper half of the complex plane, we enclose 
the very same poles whether we take iO or 2i0. The re- 
sults of the frequency integrations are thus unaltered by 
the (partial) summation of subgroup diagrams as per- 
formed above. 

Finally, we note that we have demonstrated all key 
ideas required for application to other problems. For ex- 
ample, calculations of other stationary transport quanti- 
ties (e.g. noise, adiabatic time-dependent transport) in- 
volve the same type of Kcldysh diagram a 76 i 77 and may be 
simplified exploiting the above technique. Higher order 
perturbation calculations may also come within reach. 
Importantly, the relative computational gain allowed by 
the simplifications reported here increases with the or- 



der of perturbation theory. Consider, for example, sixth 
order in Ht (neglecting the technicalities of excluding 
secular cases for simplicity). With 6 vertices to be dis- 
tributed over the two parts of the contour, taking into ac- 
count all possible contractions and directions of fermion 
lines, there are 7680 diagrams. The equivalent to Fig. [5] 
thus contains 320 irreducible and 160 reducible diagrams. 
However, one can identify 15 diagram classes comprising 
groups with x £ {0,1,2,3} vertices on the upper con- 
tour, containing 1, 1+5, 5+10, 10 diagrams respectively. 
The x = 1 groups split into one stand-alone diagram (the 
gain-loss partner of the x = diagram) plus a subgroup 
of five diagrams. The x = 2 groups comprise in turn the 
subgroup of (gain-loss partners of those) five diagrams 
plus a subgroup of ten diagrams (gain-loss partners of 
the x — 3 diagrams). Summing the diagrams in each sub- 
group and using that only one subgroup from each pair 
of gain-loss partners needs to be evaluated (cf. Sec. IV Bj) . 
our grouping method reduces the number of expressions 
to be calculated by a factor of 8. 
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